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ROOTS OF POLYNOMIALS BY RATIO OF SUCCESSIVE DERIVATIVES 
by James E. Crouse and Charles W. Putt 
Lewis Research Center 


SUMMARY 

Computer programs for finding roots of polynomials often give unsatisfactory an- 
swers where roots relatively close together are encountered. This difficulty to a large 
extent can be avoided with a procedure utilizing ratios of successive derivatives. The 
specific information gained from the ratio of successive derivatives is the number of 
roots at the root point approached and the approximate location of a trail point with re- 
spect to the closest root. The location approximation improves as a root is approached 
so a powerful convergence procedure becomes available. 

Equations are developed in this report for the general case of polynomials with com- 
plex number coefficients. Stepwise procedures are given for obtaining accurate roots 
for the general case. These principles are developed into a computer program which 
finds the real and complex number roots of polynomials for the special case of real 
number coefficients. Some examples are shown to illustrate the root resolution capabil- 
ity of the program. 


INTRODUCTION 

High-speed computing has made some formerly laborious mathematical procedures, 
such as solving for the roots of rather high degree polynomials, somewhat more prac- 
tical for broad engineering application. Before the prominence of computers the solu- 
tion of high degree polynomials for roots had an element of art to complement the sci- 
ence. The solution of general polynomials on computers, however, requires completely 
logical steps. Many methods have been developed and programmed for general use 
(e. g. , refs. 1 and 2). Almost all of these programs use iterative procedures and re- 
quire the evaluation of the polynomial at each trial root. Most of the programs work 
well for the vast majority of cases; however, they usually either compute an inaccurate 
solution or fail to converge to a solution for some root combinations. These difficulties 
are usually caused by multiple roots at a point or by two or more very close roots. 



Since roots are defined as the points for which a polynomial equals zero, iterative 
root finding techniques search for points that give a polynomial value of zero. When 
roots are not close together a polynomial will have significant slope at a root; so a tol- 
erance of the closeness of the polynomial to zero can be and effectively is used as a 
root criterion. However, when roots are very close together or when multiple roots 
occur at a point, the polynomial approaches these roots at very nearly zero slope. With 
these low slopes, very poor root resolution capability is possible with a polynomial- 
equal-zero tolerance criterion. And in some cases, programs fail to converge at all. 

In most cases of failure, either the polynomial cannot be evaluated near a root with suf- 
ficient accuracy or the polynomial value over large domains of the complex plane is 
outside the range of numbers representable on a computer. Thus, it appears that knowl- 
edge about the relative closeness of the approached root to other roots is needed for 
more comprehensive root finding computer programs. 

Polynomial derivatives give a clue as to the nature of the root or roots approached. 
In fact, the ratios of successive polynomial derivatives give the following very useful 
information: (1) the multiplicity of a root; that is, the number of roots at the root point 
approached, (2) the closeness of a trial point to the root approached, and (3) a good ap- 
proximation as to where the next nearest root is when at a root. Thus, an approach 
using the ratios of successive polynomial derivatives offers the possibility of accurate 
roots-of-polynomial computer programs with very high reliability. 

In this report the general principles of what can be learned from the ratios of poly- 
nomial derivatives (including polynomials with complex number constants) is presented 
and discussed. Then these principles are used in a computer program which finds the 
roots of polynomials for the special case of real number coefficients. The program is 
included in the report along with examples of input, output, and resolution capabilities. 


DEVELOPMENT OF THE GENERAL METHOD 
General Procedure 


In the following development of equations, the polynomial is assumed to be of the 
general form. (Symbols are defined in appendix A. ) 

2 3 i n 

P(z) = aQ + ajZ + a 2 z + a^z + . . . + a^z J + . . . a n z n (1) 


where the variable z and the constants aj for j = 0 to n are complex numbers in the 
general case. The polynomial in terms of its roots can be written as 

P(z) = (z - bj)(z - b 2 )(z - b 3 ) . . . (z - bj) . . . (z - b n ) (2) 
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Equations (1) and (2) are always analytic. Thus, a derivative always exists, and it has 
the same value at a point independent of the direction of approach. 

The first derivative of equation (2) is 

P’(z) = {[(z - b 2 )(z - bg) . . . (z - b n )] + [(z - bj)(z - bg) . . . (z - b n )] 

+ . . . + [(z -bjKz -b 2 ) . . . (z -b^j)]} (3) 

The first derivative ratio is formed by dividing equation (3) by equation (2); the re- 
sult is 


= + + . . . + _i_ (4 ) 

P(z) z - b^ z - b 2 z - bg z - b n 

Equation (4) has n terms, but note that as a root is approached a few or usually only 
one term will predominate. If only one term predominates, equation (4) can be approxi- 
mated by 


Zliil = f « _i — (5) 

P(z) z - bj 

at a trial z in the vicinity of bj. Equation (5) then can be solved for bj to give a much 
closer z approximation to the root in the next iteration. As z gets closer to bj, the 
predominance of the one term in equation (4) becomes increasingly outstanding. Thus, 
it is possible to close in on the root point very rapidly. 

The value of knowing z - bj near a root has been indicated; but so far, only the 
special case of a point near a single root has been covered. For the general case of 
arbitrary z in P(z), z may not be relatively close to any one root point so no one 
term in equation (4) predominates. Another less common possibility is the nearest root 
point may have multiple roots m, so that m terms in equation (4) sum to m/(z - bj). 
Thus, for general analysis, f ^ can be expressed as 


P'(z) _ f m 
P(z) * z - bj 




( 6 ) 


where bj is the closest root, m is number of roots at bj, and g^ is the sum of the 
rest of the terms, (n - m), in equation (4). Information about m and gj is needed 
before a value of bj can be calculated from equation (6). 

More information about m and g^ can be obtained from equations involving higher 
order derivative ratios. The development of these equations, however, becomes more 
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complicated as the order of the derivative increases; so the stepwise derivations through 
the fourth derivative are shown in appendix B. Particular equations developed there will 
be drawn into the text as needed. 

The second derivative ratio (eq. (BIO)) is 


P"(z) _ 

P'(z) 


m(m - 1) , 2mg l , 

(s-y 2 z - b i 




z 



+ gl 


( 7 ) 


This is another equation in m, g^, and a new variable gg, which is a summation of 
terms of the form l/(z - bj)^ for all bj values exclusive of the nearest root points. 
The introduction of the new variable indicates that more equations of this type will not 
give a set that can reasonably be solved in a direct mathematical way. Much implicit 
information, however, can be obtained from an order of magnitude study of the terms 
in the derivative ratios. 


Order of Magnitude Studies of Derivative Ratios 


Let us begin by observing what happens to the first and second derivative ratios 
(eqs. (6) and (7)) as a single root (m = 1) is approached. First note that as z — bj in 
equation (6), P’(z)/P(z) — ±co. 

Now consider equation (7). The first term in the numerator is zero since m = 1. 
As z — bj, equation (7) essentially reduces to 


P"(z) 

P’(z) 


2 mg. 


L = 2g 1 


m 


z - b. 


( 8 ) 


Thus, as z — bj for m = 1, the second derivative ratio approaches 2g^, which will be 
virtually a constant for small z changes near bj. Consequently, for m = 1, the first 
derivative ratio increases in magnitude rapidly and the second derivative ratio ap- 
proaches 2gj as z - bj. 

Now let us consider the case of m > 1. The first derivative ratio still approaches 
±oo as z — bj. But, in this case equation (7) essentially reduces to 


4 


m(m - 1) 


2 

P"(z) M < z ~ V _ m - 1 
P’(z) m z - bj 


( 9 ) 


Thus the second derivative ratio also approaches ±» as z — bj for m > 1. In fact, 
the ratio between the magnitude of the second and first derivative ratios is (m - l)/m 
for m > 1. 

At this point the third derivative ratio is lifted from the appendix (eq. (B15)) to show 
the following pattern that is developing: 


P ,M (z) 

P"(z) 


m(m - l)(m - 2) 3m < m “ ^1 3m ( g l " S 2 ) „ 3 

— i l + + + p ; 1 




(z - V' 


z - b i 


gl - 3 gl g 2 + 2g 3 


m(m - 1) , 2mg l , J2 „ 

— i + TT + g i" g 2 

(z - bj) 2 b j 


( 10 ) 


Consider the case of m 


2 . 


As z — b. 

J 


equation (10; essentially reduces to 


3m(m - l) gl 

P M, (z) w (z ~ b .i )2 _ 3g 

P"(z) m(m - 1) 

e ' b i )2 


HD 


As z— bj for m = 2, the third derivative ratio approaches 3gj. If m > 2, 
P’”(z)/P M (z) w (m - 2)/(z - bj as z — b.. 

J J 

Analysis of higher derivative ratios confirms the pattern indicated previously. The 
following generalizations can be made as z approaches bjt 

(1) The constant approached by the m + 1 derivative ratio is (m + l)gj. 

(2) P k+1 (z)/P k (z) — (m + 1 - k)/(z - bj) — 0 for 1 < k < m, so the number of 
derivative ratios that approach ±co is the number of roots m that are at bj. 

(3) For m > 1 and an integer k in the range 1 < k < m, 

[|p k+1 (z) |/|p k (z) |]/[l P k (z) | / 1 P k-1 (z) []- (m - k)/(m + 1 - k) where the absolute value 
symbols mean the magnitude of the vector sum of the real and imaginary parts when 
P(z) is complex. 
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These generalizations show how to interpret derivative ratios for m and g^ as a 
root is approached. When at a trial point near a root, however, it is not easy to tell if 
a derivative ratio is approaching a constant or infinity. The ratio of successive deriva- 
tive ratios, as partially introduced in generalization (3), is useful for this purpose. 

For a general integer k let us call this ratio ARATIO as it is in the computer program 
described in a later section of the report; then, 

P k+1 (z) 

ARATIO = P ^ ■ (12) 

P k (z) 

P k_1 (z) 

From generalizations (1) and (2) when k = m, 

P m+1 (z) 

ARATIO = — ^ ^ — — (m + l)g-.(z - b.) — 0 as z — b, 

P m (z) 

P m - 1 (z) 

From generalization (3), ARATIO — (m - k)/(m + 1 - k) for 1 < k < m; thus, when k 
is a positive integer less than m, | ARATIO | — C where the constant lies in the range 
0. 5 ss C < 1 and when k = m, | ARATIO | — 0. At a trial point in the vicinity of a root, 
the zero may not be very distinct; but any j ARATIO | value less than 0. 5 is an indica- 
tion that the ARATIO associated with the particular k is headed for zero. The par- 
ticular k can then be used as the current value of m. 

The values of m and obtained from a trial point provide the means of calculat- 
ing 2 - bj in equation (6) for the general case of several roots at a point. However, a 
more direct way is by equation (13) 


z 



1 

P m (z) 


(13) 


where m is both the order of derivative and the multiplicity of the root. An excellent 
characteristic of these calculated z - bj values is that they become increasingly accur- 
ate as a root is approached. They are a powerful aid in converging to a root and in es- 
tablishing a very accurate value for a root. The reason for better z - bj values as z 
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approaches b^ is the major terms in the ratio of derivative equations become increas- 
ing orders of magnitude larger than the terms ignored. 


Limitations in Practical Applications 

The preceding theoretical observations are useful only to the extent that they can 
be applied within the limitations encountered in practical work. For finding roots of 
polynomials the limitations are not severe; but they do exist; and they merit discussion. 
Almost all of the limitations are a result of the number of significant figures that can be 
carried for a constant or variable on the computer. 

The basic constants and initial parameter values that are input to the computer have 
a round off error in the last significant figure. As mathematical manipulations are 
made on the computer these round off errors and other process errors make the prob- 
able relative error of calculated parameters, such as, P(z) and its derivatives, larger. 
For meaningful ratio of derivative analysis it is necessary to recognize when the error 
of a computed value can be as large as the parameter itself. A relative error criterion 
can be established for the polynomial derivatives from an error analysis study. 

The number of significant figures that can be carried on a computer and the relative 
error criterion in essence establish the maximum size of a single meaningful derivative 
ratio. However, the judgment on the multiplicity of a trail root is made with ARATIO 
which has two derivative ratios. Thus, it is necessary to have two reasonably accurate 
derivative ratios. Since, as indicated in the earlier ARATIO discussion, both of these 
derivative ratios may be approaching infinity, the maximum allowable size of a deriva- 
tive ratio for the purpose of determining m is about the square root of the maximum 
size of a single meaningful derivative ratio; that is, about one-half the meaningful sig- 
nificant figures of a calculated derivative. This limit on the magnitude of a derivative 
ratio for the determination of m in essence establishes the minimum distance for which 
a computer program can resolve nearby roots rather than treat them as a multiple root. 
As indicated by equation (13), this minimum resolution distance for a nearby pair of 
roots is the inverse of the derivative ratio. 

When a pair of root points are resolved, the error of the root point bj is of the 
order of magnitude of the resolution distance. Usually the error of bj can be reduced 
several orders of magnitude by using equation (13) for one more iteration to obtain bj 
with the m established in resolution. The least improvement in accuracy is made if a 
pair of roots are the resolution distance apart. If a pair of roots are greater than the 
resolution distance apart, the order of magnitude of root location error reduction is the 
ratio of resolution distance to the distance between the root pair. If a root pair is less 
than the resolution distance apart, they are treated as a double root at the centroid of 
the root pair. 
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Whereas a pair of nearby roots can by resolved to a known accuracy, the resolution 
of clusters of nearby roots cannot be described as precisely. The approximate resolu- 
tion distance of an evenly spaced group of roots m on a circle in the complex plane is 
the number of meaningful significant figures of a computed derivative divided by m. 

For example, if a computer has 16-significant -figure capability, it maybe possible to 
retain about 14 significant figures in a polynomial derivative value of a tenth-degree 
polynomial. With four evenly spaced roots, the resolution distance would be only three 
and a fraction significant figures. The ratio of derivatives method, however, is most 
useful when closely packed clusters of roots or a multiroot point is encountered. In the 
approach to such a group of roots the polynomial appears to approach a high ordered 
zero or multiple root; so the actual value P(z) stays well below the absolute error as- 
sociated with a computed value of P(z) for a range of z. In the case of an actual multi- 
root point each of the m - 1 derivatives approach a lower order zero. Thus, it may 
not be possible to evaluate P(z) and its lower derivatives, but it works out nicely that 
the derivative ratios needed for root resolution (determination of m) are the ones that 
can be calculated accurately. In fact, the advantage of the derivative ratio method over 
other methods is that root analysis can still be done even though the polynomial and its 
lower order derivatives cannot be evaluated with sufficient accuracy. Upon near range 
approach to a cluster, individual roots can usually be resolved; but in the cases where 
resolution cannot be made the remaining group is treated as a multiroot located near 
the centroid of the group. 


Summary of the Algorithm 

The major features of the ratios of derivative method have been discussed at some 
length, hi the following stepwise procedure the ideas are summarized as they might be 
used to find roots of polynomials: 

(1) Find a trial z for which P(z) is in the vicinity of zero. The ratio of deriva- 
tives method usually works for this, but it may not be either efficient enough or reliable 
enough for a general program. It maybe advisable to use some standard form of two- 
or three-term Taylor’s series expansion of the polynomial for this phase. 

(2) When in the vicinity of a root evaluate the first derivative ratio, P'(z)/P(z) to 
determine the approximate location of the root. P'(z)/P(z) is an order of magnitude 
measure of the closeness of z to bj. 

(3) Find P"(z) and calculate P”(z)/P'(z). From equation (12), if P"(z)/P'(z) is 
greater than one-half of P'(z)/P(z), there is a multiple root or at least two roots within 
l/[P’(z)/P(z)] of each other so that they cannot be resolved yet. Continue taking deriv- 
ative ratios, P^(z)/P^”\z), until one is found which is less than one-half the next lower 
order one. The current multiplicity m of the root bj is k - 1 where P^(z)/P^"^(z) 
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is the first derivative ratio that is less than one-half of the next lower order derivative 
ratio. As a group of roots is approached, it may be possible to resolve roots that looked 
like multiple roots from a distance; consequently, m may be lowered during the z 
trials. 

(4) Adjust z with the following relation: 


z new 


z 


P m (z) 


P m " 1 (z) 


(14) 


As a root is more closely approached, this correction becomes better by orders of mag- 
nitude. If a multiroot point is approached in iterating, the values of P(z) and lower 
order derivatives approach high order zeros. Thus, it may not be possible to get values 
for them, but the higher derivatives can be evaluated for m and the z adjustment. 

(5) Locate z within a tolerance of about one-half the order of magnitude of the root 
resolution criterion. For the purpose of resolving nearby roots, two successive deriv- 
ative ratios which are free of round off or truncation errors are needed. Thus, a root 
resolution criterion of about one-half the significant figures that can be retained in a 
calculated derivative ratio should be used. If a z value is apparently closer to bj than 
the criterion range, it should be backed away until in the criterion range. 

(6) After z is in the root-resolution criterion range, calculate g^ from 


Sl = 


P m+1 (z) 


(15) 


P m (z) 


and improve the value of bj. When P m (z)/P m+ ^(z) is free of round off or truncation 
error, the error in bj can usually be reduced by several orders of magnitude with 


bj=Z- 


P m (z) 


P 1 "- 1 ^) 


(16) 


(7) Divide the roots out of the polynomial. The order of the polynomial will be re- 
duced from n to n - m. 

(8) Estimate the location of the next root. The value of g^ is very nearly 


Si 


1 + — L 


z - b^ z - bg 


+ . 


z - b 


(17) 


n-m 


9 



If another root is relatively close to the root just found, one term in equation ( 17) should 
predominate. A rough approximation of the initial can then be 

b.j « z - — (18) 

*1 

These general steps were used in a computer program for finding the roots of poly- 
nomials with real coefficients. A computer program for only real number polynomials 
rather than the general case of complex number polynomials is discussed for three rea- 
sons. First, far more real number polynomials are used in practical applications. 
Second, a computer program specifically constructed for only real number polynomials 
requires somewhat fewer computer operations and, thus, is more efficient for the bulk 
of problems. And third, discussion of a real number polynomial program may more 
fully illustrate the application of the ratio of derivatives concept since almost all of the 
logic needed for the general case, plus that specifically for real roots, is used. 


DESCRIPTION OF THE COMPUTER PROGRAM 

From both accuracy and efficiency considerations it is advantageous to structure a 
program to do as much analysis in real number algebra as possible. The reason is that 
fewer computer operations are required for real number computations than for complex 
number computations. Thus, this program makes a thorough search for all real roots 
first; so that the polynomials are often reduced in degree, and hence length, before 
complex number algebra is needed. The program can be considered to be composed of 
two parts, one with real algebra operations and the other with complex number opera- 
tions. The major features of each are discussed. 


Program Segments Coded in Real Number Algebra 

Most of the real number operations are done in the main program, ROOTS, which 
also serves as the control routine. The remaining real number operations are done in 
subroutines which are mentioned by name at appropriate places in the discussion. 

Preliminary calculations. - At first a check of the input data with the dimension 
limits is made. Then some tests for easy reduction of polynomial degree are made. If 
the highest degree coefficient is zero, the polynomial degree is lowered by one. The 
test is repeated until a nonzero coefficient is found. The same type of test also is made 
at the low degree end of the polynomial. A root value of zero is associated with each 
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zero coefficient on the low degree end. Thus, the roots equal to zero are immediately 
accounted for, and the polynomial degree is reduced without further ado. 

One other calculation of a preliminary nature is a possible gross scaling of the 
polynomial. Scaling makes the root resolution criterion a fraction of average root size 
instead of a fixed absolute value. The scaling is done by hexadecimal orders of magni- 
tude. Hexadecimal scaling is used since no additional error is introduced into the poly- 
nomial coefficients with computer multiplications and divisions by 16 on present com- 
puters. A scaling decision is made from a least- squares line fit of hexadecimal 
logarithms of the polynomial coefficients. If this line has a slope between +0. 5 and -0. 5, 
no gross scaling is done. 

Polynomial and polynomial derivative calculation procedure. - The real number 
polynomial can be written as 

P(x) = ajX n + a 2 x n_1 + . . . a R x + a n+1 (19) 

where the a's are now real numbers indexed from the high degree end of the polynomial. 
Indexing the polynomial coefficients in this manner conforms to the order they are input 
in the program. 

To qualitatively judge the nearness of a polynomial to a root, it is desirable to nor- 
malize equation (19). Thus, it could be rewritten as 



where the term in braces is normalized to one. As a root is approached, the term in 
braces approaches zero; but its value is meaningful only when above the absolute error 
for the calculated value of P(x). The absolute error of the term in braces, of course, 
grows with degree of the polynomial. For a tenth-degree polynomial the absolute error 
can be expected to be approximately two orders of magnitude greater than the absolute 
error of the same input number. If the first roots found are somewhat isolated, their 
relative error will be almost as good as that of P(x). In the process of dividing out a 
root, additional error is introduced into the remaining polynomial coefficients. Thus, 
even though the polynomial gets shorter, the last roots, in general, have greater rela- 
tive error. When the roots are all relatively isolated from each other, the root values 
are certainly accurate enough for engineering applications (see examples in appendix C) 
However, when several roots lie at a point or are clustered, the resolution is not 
as good as with isolated roots; consequently, the accuracy of such roots is not as good 
either. Unfortunately, when polynomials have multiple roots or clusters of roots, the 
relative error of P(x) usually is higher too. The reason is believed to be a function of 
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the relative size of the polynomial coefficients. For illustration, note that a polynomial 

of multiple roots can be written as (x - a) n . As n increases the center coefficients of 

the binomial series expansion of (x - a) n becomes orders of magnitude larger than the 

end coefficients. Another difficulty with binomial series coefficient polynomials is there 

are terms nearly equal in magnitude but of opposite size. Consider, for example, 

(x - 1)^ = x 5 - 5x^ + 10x^ - 10x^ + 5x - 1. Note that as x approaches 1, 10x** ap- 
o 

proaches lOx and the resulting sum will have few significant figures. Table I illus- 
trates how small the value of this polynomial is in the vicinity of its multiple root. 

For the purposes of illustration, table II shows how the finite precision of a com- 
puter makes it impossible to locate multiple roots accurately if the polynomial equals 
zero criterion is used. Even using ten significant figures in the evaluation of the poly- 
nomial, the value of x = 0. 999 will be accepted as a root. Thus, a ten- significant- 
figure calculation does not even yield a three- significant-figure root. Furthermore, 


TABLE n. - EVALUATION OF THE POLYNOMIAL 
(x - l) 5 at x = 0. 999 


TABLE I. - VALUES OF THE 
POLYNOMIAL (x - l) 5 AT 


VARIOUS x 

[p(x) = (x - l) 5 = x 5 - 5x 4 + 10x 3 
- 10x 2 + 5x - l] 


X 

P(x) 

0.8 

0.00032 

.9 

.00001 

.99 

. 000 000 449 

.999 

. 000 000 000 000 001 


(a) Using infinite precision 


n 

A n 

x n 

A n x n 

0 

-1 

1.0 

-1.0 

1 

5 

.999 

4. 995 

2 

-10 

. 998 001 

-9. 980 01 

3 

10 

. 997 002 999 

9.970 029 99 

4 

“5 

.996 005 996 001 

-4.980 029 980 005 

5 

1 

. 995 009 990 004 999 

. 995 009 990 004 999 


5 

P(x) = ^ A n x n = -0. 000 000 000 000 001 
n=0 


(b) Truncating to ten significant figures 


n 

A n 

x n 

v n 

0 

-1 

1.0 

-1.0 

1 

5 

. 999 

4. 995 

2 

-10 

. 998 001 

-9. 980 01 

3 

10 

.997 002 999 

9. 970 029 99 

4 

-5 

.996 005 996 0 

-4. 980 029 980 

5 

1 

. 995 009 999 00 

.995 009 990 0 



5 




P(x) = V A x n - 0. 000 000 000 0 



n=0 



12 





after dividing out the root, the coefficients of the reduced polynomial will have at most 
three significant figures. 

Again for illustration the following table shows the buildup of error as each new root 
is found with less accuracy than the preceding root: 


x 1 = 1.0063 
x 2 = 0. 9985 + 0. 0061 i 
x 3 = 0.9985 - 0. 0061 i 
x 4 = -4.9015 - 0.7623 i 
x g = -9. 4521 


5 5 2 3 2 

The polynomial, (x - 1) = x - 5x + lOx - lOx + 5x - 1, was input to a computer pro- 
gram using Laguerre's method for finding roots. The program utilizing eight significant 
figures in its calculation found the first root to about three significant figures. The next 
two roots were also located to about three significant figures but were reported to be 
complex instead of real. The last two roots are not even accurate to within one signifi- 
cant figure. However, the accuracy and internal checks within the program were satis- 
fied and the program reported no error message indicating the roots were bad. The only 
indication of trouble is the appearance of a complex root without a conjugate. 

Even with the derivative ratio program, it is better to use the most accurate calcu- 
lation procedure known for the calculation of P(x) to retain as much resolution capability 
as possible. 

An alternate method of evaluating P(x) is the following form: 



The epsilon, delta error analysis method of reference 3 applied to the previous forms 
for evaluating the polynomial P(x) shows that the latter form is probably the more ac- 
curate when clusters of roots are encountered, even though there are more operations 
used. Part of the reason is that no additional absolute error is introduced in the addi- 
tions of the exact constant one. For a fifth-degree root the relative errors of P(x) com- 
putations by equations (20) and (21) are the same; but for a tenth-degree root the relative 
error of P(x) computed by equation (21) is about 5 percent less than that of equation (20). 
There is some question as to whether this is enough difference to warrant the use of the 
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latter form; but it was used in the program since the overall program is generally quite 
fast and efficient. 

Equation (21) cannot be used directly, however, because a zero coefficient gives the 
division-by-zero problem. This problem can be handled, in general, by checking the 
magnitude of each coefficient with some tolerance which should be about the absolute 
error of other polynomial coefficients. If some polynomial coefficients a^'s are below 
the tolerance, the principle illustrated by equation (22) can be used to evaluate P(x) in 
the following way: 




Derivatives of P(x) also are polynomials, so the preceding equation forms are ap- 
plicable for them too. Calculation of the derivatives is simply a matter of setting up the 
derivative coefficients and adjusting the degree of the polynomial. The calculation of a 
polynomial and its derivatives, with x the only independent variable, is the specific 
function of subroutine RPOLY. At any point x during iteration P(x), P'(x), and P"(x) 
always are calculated. Higher order derivatives are calculated only if the ratio of de- 
rivatives segment of the program has a need for them. 

Approximate root location with second order Taylor’s series. - A second order 
Taylor’s series is used for moving x in the vicinity of a root. This method was chosen 
because, first, it is quite efficient in moving toward roots and, second, it can be pro- 
grammed to almost certainly get sufficiently near all the real roots. 

The series may be written as 

h 2 

P(x + h) = P(x) = hP’(x) + — P"(x) (23) 

2 


Since real roots are the x values for which P(x) is zero, the x increment, or h, 
which makes P(x + h) equal zero is sought. With h the only unknown in equation (23), 
the quadratic equation can be solved for values of h. If the h’s are real, the polynomial 
appears headed for a root. The value of x then is corrected by the h with the smaller 
magnitude. And P(x) and its derivatives are recalculated with the new x. When the 
absolute value of P(x) gets below the resolution distance TOLR, the ratio of derivatives 
procedure is entered. 

If the h's in the Taylor's series quadratic equation are not real, the polynomial 
either has made an approach and retreat from the x-axis or appears to be making an 
approach and retreat from the x-axis. Through logical use of P(x), P'(x) and their 
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comparison with P(x) and P*(x) from the previous x trial, the polynomial can be in- 
vestigated in this region to readily determine if this is indeed an x-axis approach and 
retreat rather than a root. If the region is only an approach, increasing x-increment 
steps are taken away from the region in the direction started. These steps are continued 
until either another real h is found or 15 steps have been taken in which case x values 
are investigated in the other direction from the starting point. If no real roots are found 
in that direction either, the complex number analysis subroutine, COMPLX, is called. 

The starting point in the search for the first real root is x = 0. The first x move- 
ment is in the direction of decreasing |P(x) | . After a root has been found, the starting 
point in the search for the next root is the estimated x value from the ratio of deriva- 
tives analysis. 

Real number analysis by ratio of derivatives . - The most important parameter in 
the program for making decisions with ratio of derivatives analysis is ARATIO, which 
is defined as 


ARATIO = 


P k+1 (x) 


P k (x) 


P k (x) 


P k_1 (x) 


(24) 


where k is the order of a polynomial derivative. This parameter is useful because we 
know from theory that 


ARATIO « 0. 5 for k = m 


(25) 


and 


1. 0 > ARATIO > 0. 5 for 1 < k < m (26) 

as a root is approached, hi the program 0. 45 was used for the 0. 5 in equations (25) 
and (26) to give a little allowance for error. These constants in essence establish when 
roots are treated as multiple or close together. 

The remaining possibility of ARATIO greater than one is an indication of either an 
x value between or among roots but not relatively near any or an encounter with absolute 
error of a computed derivative rather than a true value. To do meaningful analysis, it 
is necessary to know which it is. So it is important to establish a good polynomial error 
criterion. The criterion parameter in the program is TOLRQ. It is the probable abso- 
lute error of a computed polynomial value normalized by the lowest degree constant 
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which is nonzero. The normalized polynomial or polynomial derivative in the program 
is P(ND); and it is to be distinguished from POLY(ND), the value of the polynomial or 
polynomial derivative. Also note that since a zero subscript cannot be used in some 
computer languages, P(l) denotes the normalized polynomial, P(2) is the normalized 
first derivative, and so forth; that is, ND = k + 1. 

TOLRQ should be greater than the probable absolute error of P(ND); but excess 
margin cuts into the root resolution distance criterion TOLR, which is the square root 
of TOLRQ. Consequently, TOLRQ should be set with care. In this program TOLRQ 
varies with the degree of the polynomial as a result of an epsilon, delta error analysis 
of multiroot polynomials or varying degree. After the functional relation was estab- 
lished, the magnitude of TOLRQ was finally set by noting when computed values became 
different from known polynomial values for specific cases. With this functional rela- 
tion of TOLRQ, better root accuracy can be expected with low degree polynomials. 

The resolution criterion coefficients are essentially an absolute distance. However, 
it is probably desirable to have the resolution criterion be a relative distance to the 
local root. This effectively can be accomplished by scaling the polynomial so the local 
root is of order one. To accomplish this, subroutine SCALER is called for an order of 
magnitude check of x on only the first pass through the ratio of derivatives analysis of 
each root. If the hexadecimal logarithm of x is not within ±0.5 of zero, the polynomial 
will be hexadecimally scaled to bring the local x to order one for root analysis. How- 
ever, since the polynomial coefficients must be held within the exponent limits of a 
computer, scaling in some cases is limited. 

In recognition of the probable magnitude of polynomial absolute error, an attempt 
is made to begin the ratio of derivatives analysis with P(ND - 1), P(ND), and 
P(ND + 1) values which have magnitudes greater than TOLRQ. Initially, ND is two so 
P(l), P(2), and P(3) are calculated. The procedure is to then check P(l) with TOLRQ. 
If |P(1)| < TOLRQ, ND is increased by one if ND is less than the degree of the poly- 
nomial. Then | P(ND - 1) | is checked with TOLRQ. The procedure is repeated until 
| P(ND - 1) | is greater than TOLRQ or ND equals the degree of the polynomial. If 
|p(ND - 1) | > TOLRQ, the logical parameter LIMIT is set to zero in the program to 
indicate complete ratio of derivative analysis is possible. 

ARATIO is computed from the three polynomial and/or polynomial derivative 
values, POLY(ND - 1), POLY(ND), and POLY(ND + 1). When LIMIT = 0 and 
ARATIO < 1, analysis is done according to theory. When LIMIT = 1, the computed 
value of POLY(ND - 1) is expected to have an absolute error value rather than an ap- 
propriate value for POLY(ND - 1). Thus, the analysis methods based on an accurate 
value of POLY(ND - 1) are not used. However, some of the ratio of derivatives judg- 
ments can still be made on the assumption that the actual value of | POLY(ND - 1) | is 
probably less than the computed value of | POLY(ND - 1) | since the lowest allowable 
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program limit of |P(ND — 1) | is 0. 5 * TOLRQ. 

When ARATIO > 1 and LIMIT = 0, the present x value is between or among roots. 

Some knowledge of the nature of ARATIO contours about roots is helpful in establishing 

procedures when ARATIO is greater than one. As an example, the ARATIO contours 

o 

for the polynomial, P(z) = z (z - l)(z - 100), are shown on the complex plane in figures 
1 to 3. Such plots for known problems yield a great deal of information as to how to 
handle general cases. In this example, note that, for the group of three real roots, 

(1) the relative symmetry, (2) the points where ARATIO approaches zero and infinity, 
and (3) the values of ARATIO at derivative ratios above and below the point where an 
ARATIO approaches zero or infinity. From these figures for this single problem the 
following elements of procedure evolved: 

(1) When ARATIO for ND is greater than one, step in a direction until an ARATIO 
very near zero is found. 

(2) This low ARATIO may be a false root or roots, so check the root candidate as 
follows: 

(a) If ND = 2, either P(x) or P"(x) approaches zero; so if |P(x)| < |P M (x) |, 
x is the root. 

(b) If ND > 2, ARATIO for ND - 1 at a root should be 0. 5 by the theory. 



Figure 1. - Third degree ARATIO contours for polynomial P(z) = z 2 (z - l)(z - 100) where 
|P""(z)| /lP"’(z)l 


ARATIO = J 


IP'"(z)l / IP"(z)l 
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(3) If the root point is false (i. e. , ARATIO for ND - 1 approaches infinity when 
ND > 2 or | P"(x) | < |p(x) | when ND = 2), set x at the opposite point of symmetry and 
make the root test again. 

A summary of the basic logical steps for the ratio of derivatives analysis for real roots 
is given in appendix D. The general philosophy of resolving roots in a group is to start 
with an ND equal to m + 1 and work down in ND as root resolution progresses. The 
reason is that multiple roots should be found first. The multiple root analysis is done 
with higher derivatives, which often can be calculated accurately, whereas the lower 
derivatives and P(l) may be below TOLRQ, so that the single root analysis cannot be 
done accurately. 

If upon entry of the ratio of derivatives part of the program the x point is outside 
the group by a distance of approximately one or two times the maximum distance between 
roots in the group, the group will appear at first to be a multiple root; so the initial ND 
value naturally will be raised to the number of roots in the group plus one (note that the 
values of ARATIO away from roots in figs. 2 and 3 are always greater than 0. 5). How- 
ever, if the first x trial in the ratios of derivative analysis is within the group, it is 
probable that a low ARATIO of a false root will be found at a lower than desired ND. 
The test for this situation is to raise ND by one and check the ARATIO at the higher 
ND. If this ARATIO > 1, analysis is begun at this higher ND; but if ARATIO < 1, 
analysis is begun at the original ND. 

The aforementioned procedures probably would work quite well for only real roots; 
but a program must be able to handle the cases of real and complex roots in a group. 
ARATIO contours for known problems of this type show definite patterns; but the con- 
tours are not as distinctive as those for only real roots. These contours, however, give 
clues as to how to search for a real root when one is known to exist in the group. 

An odd number of real roots is known to exist within a range of x if P(l) is known 
to change sign in that range. As an aid for the analysis of a group of roots a running 
record of a possible root range is kept. If values of P(l) distinguishable from the abso- 
lute error are known to change sign between x trials, the logical parameter IMSURE 
is set to one. At each new x trial throughout the Taylor’s series investigation and ratio 
of derivative analysis, P(l) is checked; so the root range is often narrowed with each 
new x trial. If IMSURE is one, the root group ARATIO’ s will be quite carefully 
searched for a real root over the known root range. It may not be possible to find a root 
because the necessary P(ND) values are below TOLRQ. In this case x values are 
found which place | P( 1) | values at the ends of the root range between 0. 5 and 1. 0 times 
TOLRQ; and a root is divided out at midrange x. 

Dividing out roots . - After a root or multiple roots are identified, the normal pro- 
cedure is to (1) locate x within the range of 0. 2 * TOLR and 1. 0 * TOLR of bj, 

(2) estimate the location of the next root, (3) improve the value of the present root, and 
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(4) divide the root out of the polynomial. The x estimate for the next root is given by 
the following equation: 


x next 


x 


m + 1 



(27) 


If |P(ND - 1) | > TOLRQ, the value of the present root can be improved, usually by sev- 
eral orders of magnitudes, by the following equation: 


r r P m (x) 1 POLY(ND + 1) ] 

F^J L P0LY(ND) J 

The number of roots m found are divided out of the polynomial one at a time by the 
method presented on pages 76 and 77 of reference 3. This method is shown in appen- 
dix B for convenience. The normalized polynomial remainder (polynomial relative 
error) left after each root is divided out is also saved for output. The polynomial re- 
mainder should be zero; so its value is a good clue as to whether or not the root was 
determined as accurately as it should have been. 

When the degree of the polynomial is reduced to two or less, subroutine QED is 
called. Any remaining roots are found by direct computation. Then all of the polyno- 
mial roots and remainders are printed. 


Program Segments in Coded in Complex Number Algebra 

If more than two roots remain after the search for real roots has been completed, 
subroutine COMPLX is called for the search of complex conjugate pairs of roots. The 
procedures in COMPLX somewhat parallels those of ROOTS. 

Polynomial and polynomial derivative calculation procedure. - The polynomial P(z) 
and its derivatives are calculated with a form of equation (21) for complex numbers. The 
calculation of a polynomial and its derivatives with z the independent variable is the 
specific function of subroutine CPOLY. At any z point during iteration, P(z), P’(z), 
and P"(z) always are calculated. Higher order derivatives are calculated only if the 
ratio of derivatives segment of the program has a need for them. 

Approximate root location with second order Taylor’s series. - A second order 
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complex Taylor’s series is used to locate z in the vicinity of a root. The series may 
be written as 


P(z r ) = P(z) + P’(z)(z r - z) + P”(z) 


(z„ - z)' 


(29) 


Since roots are z r values from which P(z r ) is zero, equation (29) can be written as 


0 = P(x) + iP(y) + jP’(x) + iP'(y)J(z r - z) + {P"(x) + iP"(yj] 


(z„ - z) 2 


(30) 


The solution of equation (30) for z r - z can be expressed as 


a ■ , _ -Qp’(x) + P’(y)] 

r P"(x) + iP”(y) 


s ’(x)] 2 - [P’(y)] 2 - 2[p(x)p"(x) - P(y)P"(y)] + 2i [P’(x)P’(y) - P(y)P”(x) - P(x)P”(y)] 

P"(x) + iP"(y) 


(31) 


The plus or minus sign on the square root term of equation (31) indicates that the term 
may be used as computed or at a 180° phase angle. Both solutions for z r - z are com- 
puted. The one that gives the smaller absolute value of z r - z in the complex plane is 
used for the new trial z r> When the absolute value of P(z p ) is reduced below the reso- 
lution tolerance TOLI, the ratio of derivatives analysis is begun. 

The initial coordinates of the Taylor's series search for the first complex root are 
the x that corresponds to some local minimum of |P(x)| in the real number Taylor's 
series analysis in ROOTS and a y value equal to the square root of | P(x) | . After a 
complex conjugate pair of roots have been found, the initial coordinates in the search 
for the next roots are the estimates from the final stages of the ratio of derivatives 
analysis of the present root. 

Comp lex number root analysis by ratio of derivatives . - It is rather interesting that 
all of the information extracted from real number ratio of derivative analysis along the 
x-axis can also be extracted from complex number ratio of derivative analysis in the 
complex plane. The reason for this lies with the very definition of analytic complex 
number derivatives. The essential point is that a derivative has the same value at a 
point no matter from which direction the point is approached. The complex derivatives 
have real and complex parts which can be vectorially combined into derivative magni- 
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tudes. Ratio of derivative analysis based on these derivative magnitudes (ARATIO'S) 
works in the same way as real number ratio of derivative analysis for the determination 
of m. 

As with real roots the derivative ratios used to establish m are the ones from 
which the component adjustments are made. However, in complex numbers it is neces- 
sary to break down the appropriate derivative ratios into real and imaginary parts to 
get x and y adjustment values. Since the development of appropriate equations is 
somewhat complicated, it is done in appendix B, and only the results are shown here. 
The adjusted or new x and y values are obtained from the derivatives by the following 
equations: 


x new " x 


A 2 + B 2 


(B32) 


y = v + 

J new J 


B 


A 2 + B 2 


(B33) 


where 


m-l/..\ . T>m,„\T>m-l( y ) 


A _ P m (x)P m ~ A (x) + P (y)P 

[P m - 1 (y >] 2 


(B27) 


and 


E _ P (y)P (x) - P (x)P (y) 

[p m - 1 (x)] 2 + [p m - 1 (y>] 2 


(B28) 


In the ratio of derivatives analysis, the logical parameter LIMIT is again used to 
indicate whether or not the computed polynomial value is greater than the probable ab- 
solute error. The complex number criterion TOLIQ is of the same form as TOLRQ; 
but TOLIQ was set at four times TOLRQ. This is because more computer operations 
are needed to compute a complex number polynomial value than for a real number poly- 
nomial of the same degree. The complex root resolution criterion TOLI is the square 
root of TOLIQ. 

For the initial z trial in ratio of derivatives analysis two checks are made. First, 
the degree of the root approached is raised as necessary within the restraints of the de- 
gree of the polynomial to begin the analysis with LIMIT equal to zero if possible. And 
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second, the polynomial is scaled if the magnitude of z is not sufficiently near one. 

If, during ratio of derivatives analysis, ARATIO is less than one and LIMIT is 
zero, the analysis follows the theory. When ARATIO is greater than one with LIMIT 
equal to zero, the present point lies between or among roots. Since there is another 
adjustment degree of freedom on a plane as compared to that along a line, there is 
greater desire to have z movement procedures which effectively and efficiently move 
toward roots. Again the ARATIO contours in the complex plane give clues of effective 
procedures. See figure 4 for an example of contours about complex roots. Also, the 
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Figure 4. - ARATIO contours for polynomial P(z) * (z - 9 + i) 2 (z - 9 - i) 2 (z - 11 + i) 2 (z - 11 - i) 2 . 
(Note: The polynomial is symmetric in four quadrants about the point x * 10, y = 0. ) 


ARATIO contours of figures 1 to 3 are quite similar to those in the vicinity of a similar 
grouping of complex roots. When ARATIO is somewhat greater than one, the contours 
are nearly circles. Clearly steps should be toward lower ARATIO; but a direction 
started may not necessarily head in the direction of an ARATIO less than 0. 45. Con- 
sequently, the need of a curved path is indicated. 

In the program the procedure, when ARATIO is greater than one, is to test step 
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in the x direction and then in the y direction. From these two normal steps the most 
effective direction for reducing ARATIO is determined. Then a spiral path in 30° in- 
crements is followed in the search for an ARATIO less than 0. 45. If while on the 
spiral ARATIO increases from the previous iteration, analysis reverts back to the pre- 
vious point. At that point, a new direction is determined from another set of x and y 
test steps; and then a spiral turning in the opposite direction is begun. The experience 
is that the program should seldom have to move more than seven points on any one 
spiral or should seldom have to begin a new spiral more than once. 

In the complex plane the general philosophy of resolving roots in a group is nearly 
the same as on the real axis in ROOTS. There are some differences, however, which 
deserve discussion. First, since there is freedom to move about the complex plane, 
the problem of resolving complex from real roots no longer exists. If it turns out that 
there still are real roots in the group, they can be identified in the complex plane too. 

A second difference is in the check of the degree of the root candidate investigated. 

In ROOTS a root candidate was quite well located before a final check of the root degree 
(ND level) was made. In the complex plane there are not quite as many logical param- 
eter possibilities, so it is a little easier to change the root degree at any time during 
the ratio of derivative analysis. This is done by checking ARATIO and another 
ARATIO at the next lower ND for each z trial when the degree of the root is greater 
than one. 

A third difference is that the symmetry of ARATIO contours cannot be depended 
upon. In a root group it is less probable that points where ARATIO equals zero will 
be symmetrical with each other about a point where ARATIO equals infinity. However, 
by keeping track of ARATIO* s at the two ND's it is usually possible to identify false 
roots quickly. 

Dividing out roots . - As a root or multiple roots are approached, the final decision 
on the multiplicity of the root is made when z is within the range of 0. 2 * TOLI and 
1. 0 * TOLI of bj. After m is established, the location of the next nearest root is 
estimated from gj found by 


(3 

m + i p m (z) 

Since the conjugate of the present root is also in gp the breakdown of gj for the next 
nearest root location is a rather lengthy development which is shown in appendix B. 

The resulting equations for the new coordinates are 
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( 37 ) 


and 



where g R and gp respectively, are the real and imaginary parts of g^. The values 
of the present root coordinates then are improved with the same equations used for 
coordinate adjustment, that is, equations (32) and (33). 

A complex root is first artificially divided out of the polynomial to find the real 
and imaginary normalized polynomial remainders. Since these remainders should be 
zero, their values are clues as to whether or not the complex root was located as ac- 
curately as it should have been. The dividing out process is artificial in the sense that 
the polynomial coefficients are not permanently changed until the root and its conjugate 
are divided out together in real algebra. The equations and details for dividing out the 
complex roots by both methods are shown in appendix B. The polynomial remainder 
from the real algebra division appears in the same printout line with the conjugate root. 


Examples 

The general capabilities and limitations of the program have been discussed; but 
they are most effectively shown with examples. The exact roots of several polynomials 
that are difficult to solve are shown with the computed roots in the program output sec- 
tion of appendix C. These examples illustrate the root resolution capability of this pro- 
gram. 

Appendix C also contains a listing of the program, a description list of the program 
variables, and other special instructions for a user. 


CONCLUDING REMARKS 

The ratio of derivatives method is a powerful method of finding roots of polynomials. 
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A computer program for finding the roots of real number polynomials by this method was 
developed to see how useful the theory is within the confines of calculated number accur- 
acies. Examples of the root resolution power of this program illustrate that the method 
is indeed powerful in practice as well as in theory. Although a program was not devel- 
oped for complex number polynomials coefficients, equally effective root finding pro- 
grams can be developed with the ratio of derivatives method. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 2, 1972, 

764-74. 
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APPENDIX A 


A 

a 

B 

b 

c 


*1 


Sk 

h 

k 

l 

m 

n 

P(x) 

P(z) 

P k (z) 

x 

y 

z 


SYMBOLS 

combination of complex derivative terms (see eq. (B27)) 
general polynomial coefficient 

combination of complex derivative terms (see eq. (B28)) 
general coordinate of roots 

general polynomial coefficient after a root is divided out 

n 

function of the form ^ [l/(z - b-)J 

j=l J 

function of the form ^ \}/( z ~ bj) k J 

3=1 

n-m 

f 1 at a root with the local root point excluded, fl/(z - b-)l 

j = l 3 

n-m 

f^ at a root with the local root point excluded, ^ [l/(z - bj) k J 

x adjustment increment toward a root by Taylor's series 

degree of an arbitrary polynomial derivative 

arbitrary exponent of the function fj in eq. (B21) 

number of roots at the particular root point 

number of roots in the polynomial 

polynomial in real numbers 

polynomial in complex numbers 

k in polynomial derivative in complex numbers 
independent variable or real component of independent variable 
imaginary component of independent variable 
independent variable in complex number form 


Subscripts: 

I refers to imaginary part of a complex variable 

j an arbitrary term or root of a polynomial 



■ I 


n n** 1 and last root or next to last polynomial coefficient 

new newest approximation to a root 

next estimate of next root location 

R refers to real part of a complex variable 

r root or very near root value 

Superscripts: 

j arbitrary term in the polynomial 

k degree of an arbitrary polynomial derivative 
l arbitrary power of the function f^ in eq. (B21) 

m number of roots at a root point 

n degree of the polynomial 

' first derivative 

" second derivative 

third derivative 
"" fourth derivative 
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APPENDIX B 


RATIO OF POLYNOMIAL DERIVATIVE EQUATIONS 
Development of General Equations 


A general complex number polynomial in terms of its roots can be written as 
P(z) = (z - b 1 )(z - b 2 )(z - b 3 ) . . . (z - b n ) 

The first derivative of equation (Bl) is 

P’(z) = {[(z - b 2 )(z - b 3 ) . . . (z - b n )] + [(z - bjXz - bg) . . . (z - b R )] 

+ • • • + [(z - bj)(z - b 2 ) . . . (z - b n _ 1 )]J 
Form the first derivative ratio by dividing equation (B2) by equation (Bl) 

= ! — + — 1 — + — 1 — + . . . + _1_ + . . . + _i_ 

P(z) z - b^ z - b 2 z - b 3 z - bj z - b fl 

Generalize equation (B3) to 


f 


1 " 


m 


- b i 




where bj is the closest root to z, m is the number of roots at bj, and g^ is the 
of the remaining n - m terms in equation (B3). Equation (B2) can be written as 


P’(z) = P(z) • = P(z) • f, 

P(z) 1 


Differentiate equation (B5) to get the second derivative of P(x) 


P"(z) = P’(z) • f 1 + P(z) • f’j 


= P(z) • i[ - P(z) 


(z - b x ) 2 (z - b 2 ) 2 (z - b g ) 2 


(z - V 


(Bl) 


(B2) 


(B3) 


(B4) 

sum 

(B5) 
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P"(z) = P(z) • (f 2 - f 2 ) 


(B6) 


where 


f 2 ” 


- 1 — + L_ 

_(z - b x ) 2 (z - b 


/ < z - b 3> 2 


+ •••+■ 


(z 


J. 

- b ) 2 
n' J 


Generalize equation (B7) to 


f 2 " 


m 


(Z - b,)‘ 


62 


The second derivative ratio is 


P"(z) _ f l ~ f 2 
P f (z) fj 

Substituting equations (B4) and (B8), equation (B9) becomes 


P M (z) _ 
P’(z) 



m(m - 1) [ 2mg l 

(z-b/ 2 ~ b j 


+ g 2 l~S 2 


z 


m 



+ *1 


Differentiate equation (B6) to get the third derivative of P(z) 

P"’(z) = P’(z) • (ff - f 2 ) + P(z) • (2f 1 f , 1 - f^) 

= P(z) • fj - (f 2 - f 2 ) + P(z) • [tf^-Dfg - (-2)f 3 ] 
= P(z) • (f 2 - 3f 1 f 2 + 2f 3 ) 


(B7) 


(B8) 


(B9) 


(BIO) 


(Bll) 
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where 


f 3 " 


1 + -i + — -1 


(2 - b x ) 3 (2 - bg) 3 (2 - bg) 3 


+ • . • + 


(z 


-V s . 


Generali2e equation (B12) to 


(B12) 


f o — 


m 


(z - b.r 


S3 


(B13) 


The third derivative ratio is 


P M, (z) = f l ~ 3f l f 2 + 2f 3 


f2 _ f 


P"(z) 

With the generalized f substitutions (B14) becomes 


(B14) 


P M, (z) _ 
P M (z) 


m(m - l)(m - 2) . 3m < m " . 3m Sl . „3 

+ + r g^ 


(z ~ V 


(i-V 


z - b. 


3g 2(r J2 r + g i) +2g s 


i z - b- 
^ J 


m(m ■ 1) , 2n *l , 2 


(z-bjV 


z - b. 
1 


gl -g 2 


(B15) 


Differentiate equation (Bll) to get the fourth derivative of P(z) 

P ,,M (z) = P*(z) • (f| - 3f jf 2 + 2f s ) + P(z) • (3^ - 3f 2 f , 1 - 3f x f 2 + 2fj) 

= p(z). t v (f 3 -3f 1 f 2 + 2fg) +P(z).|3(f 3 - f 2 )(-l)f 2 - 3f — 2)f 3 + 2(-3)f 4 

= P(z) • (4 - 6f}f 2 + 8f jfg + 3fJ - 6f 4 ) 

where 


(B16) 


f 4 = 


(z - b x ) 4 (z - b 2 ) 4 (z - bg ) 4 


+ . . • + ' 


(z - b n ) 


(B17) 
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Again generalize the f as before 


f 4 = 


m 


(z - V 


g 4 


(B18) 


The fourth derivative ratio is 


p,',’(z) 4 - 6f l f 2 * 8f l f 3 * 3f 2 - 6t 4 


q - 3^2 + 2f 3 


P M, (z) 

With the generalized f substitutions, equation (B 19) becomes 


(B19) 


m(m - l)(m - 2)(m - 3) 4m(m " ~ 2 ^1 6m ( m “ 4m Sl 4 

1 Ci L + + + + g^ 


P" M (z) _ L ( z ~ b j) 

P"’(z) 


f - V 


(z-bj 2 z - b J 


mtm - l)(m - 2) + 3 m ( m ~ + ^ m ^l + ^3 




(z - bj)' 


Z-bj 


3 gc 


m 


v z ' b j 


+ gl +2g 3 


“6gf 


m(m - 1) , 2m Sl r g 2 


(z - b-) 2 2 - b j 


+ 8g3 (r^ +gi ) +3g 2- 6g 4 


m(m - l)(m - 2) , 3 m < m ~ ^gl , 3 mg l , g 3 


L (* - b / 


( Z - b/ 


z-^ 


■ 3B2 (r^ + Sl ) + 2g3 

(B20) 


Higher derivatives develop in the same way. The process follows the pattern used to get 
the first four derivatives. The higher derivatives have progressively more terms in the 
generalized forms for the f's. This does not seem consistent with the fact that P n (z) 
must be a constant and all the higher derivatives must be zero. But what is happening is 
that in progressing to high derivatives with the generalized forms for the f's, more and 
more terms are carried internally that would cancel if the f’s and g’s were expanded in 
their z - b terms. 

At this point a generalization of the patterns shown by the first four derivatives is 
in order. First note that the terms in the brackets in equation (B20) and in each of the 
corresponding lower order equations correspond to something like a binomial expansion. 

J.U 

The general form for fj to the k tn power can be expressed as 
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(B21) 


k 



k! m! (g x ) Z 

i!(k - Z)!(m - k + l)l(z - b-) k_i 


if O', is understood to be one. 

Comparison of equations (B19) and (B20) and each of the other corresponding pairs 
of derivative ratio equations indicates the following procedure for transferring the deriv- 
ative ratio equations from those with f's to the form with m's and g’s: (1) replace the 
f^'s and powers of fj with equation (B21), and (2) replace all the f’s that do not have a 
subscript of 1 with g's that have the corresponding subscript. 


Division of a Real Root from the Polynomials (from ref. 3) 


The general polynomial 

P(x) = a^x 11 + a 2 x n_1 + . . . + a n x + «. n+ j 
can be expressed in terms of a root as 

P(x) = (x - x^^x 11 - 1 + c 2 x n " 2 + . . . + c n _jx + c n ) + c n+1 
By equating powers of x in equations (B22) and (B23) 


(B22) 


(B23) 


a l = c l 


a 2 = c 2 - c l x r 


a 3 “ c 3 _ c 2 x r 


a = c - c i x 
n n n-1 r 


a . = c - c x 
n+1 n+1 n r 


(B24) 


The c's can be solved for directly by starting from the top of the equation set (B24) 
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c l = a l 


e 2 = a 2 + c lXr 
e 3 = a 3 + c 2 x r 


c n+l = a n+l + c n x r 

The constants, Cj to c n , are the coefficients for the polynomial left after the root is 
divided out. The constant c n+ ^ is the remainder which should be zero. 


Complex Number Root Coordinate Adjustment by Ratio of Derivatives 

When in the vicinity of roots, complex derivative ratios are used for the analysis 
and finer adjustment to root coordinates. Once m is determined, the coordinate adjust- 
ment equations are obtained from the following general equations: 

■ P & * (B26) 

P m-1 (z) z " z new 

Break equation (B26) into real and imaginary parts 


[p m (x) + iP m (yj] Tp™ - 1 ^) - iP m ~ 1 (y)] _ 1 P x x new^ ~ y newO 

[p“-l(x) + iP m ~*(y)] - iP m_ ^(y)J C x ~ x new^ + “ y newO C x ” x new^ " ~ ^newO 


P m (x)P m_ 1 (x) + P m (y)P m ” 1 (y) + i _P m (y)P m l ( x ) ~ P m (x)P m 1 (y)j = ^ x x nevP y new^ 


[p 111 ' 1 ^)] + jp m " 1 (y)J 


< x - x new )2 + ^ - y new )2 


Let 


A P m (x)P m ~ 1 (x) + P m (y)P m ' 1 (y) 

[p m_ 1 (x)] 2 + [p m_ 1 (y )] 2 


(B27) 
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(B28) 


Therefore, 


B _ P m (y)P m ~ 1 (x) - P m (x)P m " 1 (y) 
[p m_1 (x)] 2 + [P m_1 (y)] 2 


A + iB = 


< x - x new> - «y - W 

“ x new^ + (y “ ynew^ 


A = 


X - X. 


new 


( x “ x new^ + ^ ” ynew^ 


(B29) 


B = - 


y - y, 


new 


^ x " x new^ + ^ ynew^ 


(B30) 


By dividing equation (B29) by equation (B30), we have the following simple relation: 


A 

B 


x - 


y - 


x new 

V 

J new 


(B31) 


Equation (B31) then can be used in equations (B29) and (B30) to get the following x and 
y adjustment equations: 


x new = x ' 


A 2 + B 2 


y = y + — — 

J new J o o 

A 2 + B 2 


Estimate of Next Nearest Complex Root 

When near a root, the term g^ is obtained from 

„ _ 1 P m+1 (z) 

g l , 

m + 1 P (z) 


(B32) 


(B33) 
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where 


n-m 


*1 


=7 — 

ti z - b i 


Included in the summation is the conjugate of the present root. If the next nearest root 
and its conjugate are relatively much nearer than the others remaining, g ^ can be ap- 
proximated by 


. m 1 , l 

gl = g R + ig x = =- + + = 

z-z z-z , z - z . 
r r r next r next 

where g R and gj are, respectively, the real and imaginary parts of g^. Now continue 
to expand the preceding equation as follows:, 


g R + igj = 


m 


(x r + iy r> - (x r " iy r> < x r + ^ ' ( x next ' iy next> < x r + ly r> " (x next " iy next> 


2iy r + (x r - x next ) + i(y r - y next ) + (x r - x next ) + i(y r + y next ) 

-im , [ (x r ~ x next> + i(y r + ^exp] + [ (x r ~ x next> + *< y r ' y next>] 

2y r [(* r - x next ) + i(y r - y next )] [(x r - x next ) + i(y r + y next )] 

^im + 2 < x r - x next> + 2iy r 

2Yr < x r - x next> 2 ' ( y r " y next) + 2i < x r ' x next )y r 

-im + 2 B x r x nexP + ^ y 3 [^ x r ~ x nexP ~ ~ y next) ~ 2 ‘^ x r ~ x nexP y rJ 

2yr [ (x r - x next )2 “ ( y r ' y next)] 2 + 4(x r ‘ x next> 2y r 

-im , 2 fr x r ~ x next>[( x r ~ x next> 2 ~ ( y r ~ y next)] + 2 < x r ~ x next )y r} 

2y r D 

, 2iy rf[ (x r - x next )2 ~ ( y r ~ y next)] ~ 2(x r ~ x nex/) 
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where 


D = f x r " x next> 2 " (4 ~ ^ext)] 2 + ^ x r ~ x nexA 


The separate equations for real and imaginary parts are 


g R ” 


g I 


2 < x r - x next>fr x r ' x nex/ + (y 2 + 4eJ\ 


D 


+ m _ _ tyr[f x r ~ x nexP + (yr ~ ynext)] 


2y T 


D 


Divide the preceding equations and solve for y nex ^ 


g I 


m 

2y. 


< x r “ x next> 

[( x r " x next 

) 2 + | 

{4 + y 2 ext)] 

y r 

h- 

' x next )2 + 1 

(4- 

y 2 Y 

J next/. 



' g R y r [< x r - x nex/ + (*? " y^ext)] = + ^ x r ' x nexP [< x r “ x nex/ + 4 + 5 


v 2 

J next 


g Ry r " ( g I + — j^ x r - x next> 


2 

ynext 


|g R y r + + 


* [<*r - x next > 2 + *?] + " x next 

>] [< x r - *next )2 + *?] 


— - * 


r next' 


g Ry r ‘ ( g I + ^ K x p - x next> 


Substitute equation (B36) into equation (B34) and solve for x neX j- 


(B34) 


(B35) 


nextj 


2 *y| 


(B36) 
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I 


x next = x r 


gRy r 

4 + fc + — 
l 1 


(B37) 


+ St + 


m 

2y„ 


The plus y next value is used in equation (B36) in the attempt to keep the search for 
complex roots above the x-axis. 


Division of a Complex Root from the Polynomial 


The general polynomial 

P(z) = a^z 11 + a 2 z n_1 + . . . a n z + a R+ ^ (B38) 

can be expressed in terms of a root as 

P(z) = (z - z r ) (cjz 11-1 + c 2 z n “ 2 + . . . c n _jz + c n ) + c n+1 
Equating powers of z results in the following: 


a l " C 1 

a 2 = c 2 " c l z r 
a 3 = c 3 - c 2 z r 

f ( B3£ >) 


2, — C — C i Z 

n n n-1 r 

a i=c - c z 
n+l n+1 nr J 

The c's are complex coefficients which may be solved for directly by starting from the 
top of the equation set (B39). 
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C 1 - a l 


c 2 = a 2 + c l z r 
c 3 = a 3 + c 2 z r 


c n+l a n+l + c n z r 
so 

C R, 1 = a l 

C R, 2 = a 2 + C R, l x r ' C I, l y r 
C R, 3 = a 3 + C R, 2 x r " C I, 2 y r 

C R, n+1 = a n+l + C R, n x r ' C I, n y r 

and 

C I, 1 = 0 

C I, 2 " C R, l y r + C I, l x r 
C I, 3 - C R, 2 y r + C I, 2 x r 


1, n+1 


c C R, n^r + C I, n x r 


The sets of constants c R i to c r n and c j i to c j n ) are the respective real and 
imaginary sets of coefficients for the polynomial left after the root is divided out. This 


division is used only to determine the remainder terms c^ n+ ^ and 


T, n+1' 
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Division of a Complex Conjugate Pair of Roots from the Polynomial 

The general polynomial (eq. (B38)) can be expressed in terms of a root and its com- 
plement as 

P(z) = (z - z r )(z - z r ) (c^z 11-2 + c 2 z n ' 3 + . . . + c n _ 2 z + c n _j) + c n z + c n+1 

The product, (z - z r )(z - z r ), can be expressed as 

(z - z r )(z - z r ) = z 2 - z(z r + z r ) + z r z r 

= (x + iy) 2 - (x + iy)2x r + (x r + iy r )(x r - iy r ) 

= x 2 - y 2 - 2xx r + x 2 + y 2 + 2iy(x - x r ) (B41) 

Since the problem being solved is by definition a real polynomial with x the only inde- 
pendent variable, a general y does not exist. Therefore, equation (B41) can be written 
in real algebra. 

(z - z r )(z - z r ) = x 2 - 2xx r + x 2 + y 2 
The general polynomial can be written as 

P(x) = (x 2 - 2xx r + x 2 + y 2 ) (cjx 11-2 + c 2 x n “ 3 + ...c n .j)+ c n x + c n+1 

Equating powers of x 
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a l - C 1 

a 2 = e 2 - 2 Cl x r 

a 3 - c 3 - 2c 2 x r + c l( x r + y ?) 
a 4 * c 4 - 2c 3 x r + c 2( x r + y ?) 


(B42) 


a n ' c n - 2c n-l x r + c n-2( x ? + y r) 
a n+l - c n+l + c n-l( x r + y ?) 

The equation set (B42) can be solved directly for the c’s by starting from the top. 

C 1 = a l 

c 2 = a 2 + 2c i x r 

c 3 - a 3 + 2c 2 x r - c l( x i- + y ?) 

c 4 = a 4 + 2c 3 x r - c 2 (x 2 + yj) (B43) 

c n - a n + 2c n-l x r ' c n-2( x ? + y r) 
c n+l = Vl - c n-l( x r + y r) 

In this case the polynomial remainder after the root and its complement are divided out 
is 

c n x r + c n+l (B44) 
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APPENDIX C 


THE PROGRAM 
Information for Users 

The computer program is written in FORTRAN IV language and it is run in double 
precision. On an IBM direct coupled 7044 - 7094 system the running time averages 
about 0. 01 minute for an eighth-degree polynomial. The program has one systems sub- 
routine called DUBIO which appears before the input READ statement. The function of 
this subroutine is to allow double precision input and output formats on the NASA Lewis 
computer. It is not needed on most other computers. For some applications it may not 
be necessary to use the double precision input and output, but one should always remem- 
ber that the output accuracy is a function of input accuracy. 

The first card of an input data set has the number N of polynomial coefficients ex- 
pected on the following data cards. The polynomial coefficients A(I) are read in order 
beginning with the high degree end of the polynomial. The input format is shown in fig- 
ure 5. 



Figure 5. - Input data. 
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The polynomial relative error criterion constants, TOLRQ and TOLIQ, which are 
defined near the beginning of the main program ROOTS, are based on 16-significant- 
figure double precision capability. If the program is to be used with other than 
16 -significant -figure machine capability, the constants in TOLRQ and TOLIQ should 
be adjusted directly with the change in significant figure capability. 

Note that the way the program is structured, meaningful analysis can only be done 
if TOLRQ and TOLIQ are above the absolute error levels of computed polynomial 
values. However, raising TOLRQ and TOLIQ decreases the root resolution capability 
of a program. To allow for near maximum program root resolution capability, TOLRQ 
and TOLIQ are internally adjusted to account for the greater number of computations 
needed to evaluate a higher degree polynomial. In the Polynomial and Polynomial De- 
rivative Evaluation section’s discussion of ROOTS it was pointed out that the least capa- 
bility occurs when roots are grouped. The TOLRQ and TOLIQ adjustments with 
polynomial degree are based on roots being grouped up to degree ten. For higher degree 
polynomials the probability of groups of more than ten roots is small, so the increase of 
TOLRQ and TOLIQ with polynomial degree were somewhat leveled off. 

A final comment concerns the use of ARATIO contour plots if a user is not satis- 
fied with the computed roots or has reason to believe that a ratio of derivatives program 
is not working properly for a particular problem. The polynomial coefficients can be 
used in a little program to compute ARATIO values for two or three ND values on a 
x - y grid over the troublesome area. Contour plots from the grid values can be a very 
helpful aid in indicating where the roots should be. Values of POLY(ND - 1) in the ab- 
solute error range will fog the issue, but such ranges are usually obvious on a plot. 


Description of Program Variables 


Symbol 


Description 


A(I) 

ARATIO 

B(I) 

CA(I) 

CRATIO(I) 

D 


polynomial coefficients beginning from high degree end 

absolute value of the ratio of successive ND derivative ratios, 

[| POLY(ND + 1) | / 1 POLY(ND) |]/ [| POLY(ND) | / 1 POLY(ND - 1 ) |] 

coefficients of polynomial derivative 

polynomial normalizing constant; also the calculated imaginary part of 
polynomial coefficients at a complex root 

imaginary part of a complex number derivative ratio 

value of the complex polynomial, POLY(ND - 1), squared 
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DELTX 

DELTY 

DIR 

DRATIO 

DTHETA 

DX 

DXT 

DYT 

DZT 

DZT2 

D1 

D2 

D3 

FACT 

FK 

H 

HX 

HY 

HYP 

HI 

I 

ICHECK 


real component of an independent variable step in the complex plane 

imaginary component of an independent variable step in the complex plane 

step direction indicator in the Taylor's series search for real roots 

ratio of successive ND derivative ratios, [POLY(ND + 1) /POLY(ND)]^ 
[pOLY(ND)/POLY(ND - 1)] ' 

spiral angular increment (30°) used in moving toward a root when known to 
be between nearby roots in the complex plane 


x-increment to estimated location of next root; also, a reference x step 
when searching for a root known to be within a root group 

normalized real part of complex number polynomial first derivative as 
used in complex number Taylor's series solution for a new z 


normalized imaginary part of complex number polynomial first derivative 
as used in complex number Taylor’s series solution for a new z 


relative distance to root with plus sign in complex number Taylor’s series 
equation 


relative distance to root with minus sign in complex number Taylor's 
series equation 

value of complex polynomial, POLY(ND), squared 
value of the ND complex derivative ratio squared 
value of the ND - 1 complex derivative ratio squared 


ratio of x adjustment to previous x adjustment in Taylor’s series ap- 
proach to a root 

multiple of coefficient for calculation of polynomial derivative 

x adjustment toward a root by the Taylor’s series approximation 

x adjustment toward a root by complex Taylor’s series approximation 

y adjustment toward a root by complex Taylor’s series approximation 

magnitude of the square root term in complex plane Taylor's series approx- 
imation for a new z 


x backsteps in Taylor’s series search for real roots 

polynomial coefficient subscript for dimensional variables 

a logical parameter which activates a running record and update of XHIGH 
and XLOW 
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DDIR 

n 

IMSURE 

INK 

IR 

IREAL 

ISIGN 
IT RIG 

IW 

IXM 

J 

JTRIG 

K 

KTRIG 

L 

LIMIT 

LTRIG 

M 

MULT 

N 

NC 

ND 

NDEG 

NDRV 


logical direction indicator in the Taylor’s series search for real roots 

index of previous nonzero polynomial coefficient in the computation of a 
polynomial value 

logical indicator of the certainty of a real root remaining in the polynomial 
logical indicator of the direction of scaling 
system input tape number 

logical indicator of the type of root sought (IREAL = 1 for real roots, 
IREAL = 0 for complex roots) 

routing device for setting spiral direction 

routing device of how polynomial approaches the x-axis in the Taylor's 
series search for real roots 

system output tape number 

logical indicator used when a pair of root candidates need to be investigated 
current number of roots found 
routing device 

index used in dividing out roots 

routing device in ratio of derivative analysis for roots that are difficult to 
resolve 

logical device for counting the number of successive zero coefficients in a 
polynomial 

logical indicator that a computed polynomial value is probably in the abso- 
lute error range 

logical device and counter used when known to be between or among roots 
degree of polynomial derivative 
number of roots at a root point 

number of polynomial coefficients (degree of the polynomial plus one) 

index used in scaling a polynomial, also a counter during preliminary 
checks 

highest order polynomial derivative needed in current analysis 
degree of the current polynomial 
order of the polynomial derivative 
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NEXP 

NEXPS 

NK 

NM 

NND 

NUM 

OARAT 

OAX 

OAY 

ODRAT 

OPOLY1 

OPOLY2 

ORATIO 

PHI 

P(I) 

POLY(I) 

POLY1R 

PXHIGH 

PXLOW 

Q 

Ql 

RA 

RA(I) 

RAD 

RATICKI) 


exponent hexadecimal order of magnitude factor used in polynomial scaling 

exponent hexadecimal order of magnitude factor used in initial polynomial 
scaling 

index used in polynomial scaling 

index used in scaling a polynomial, also a counter during preliminary 
checks 

temporary storages of ND 
number of roots in the polynomial 
value of ARATIO for previous iteration 

reference value of ARATIO for the x-step in setting up the spiral refer- 
ence angle 

value of ARATIO from the x-step in setting up the spiral reference angle 

value of DRATIO for previous iteration 

most recent good value of POLY(l) 

most recent good value of POLY(2) 

value of RATIO(ND - 1) from previous iteration 

complex plane angle of square root term in complex number Taylor's 
series solution for a new z 

normalized polynomial, P(l), or a polynomial derivative when I > 1 

polynomial, POLY(l), or a polynomial derivative when I > 1 

ratio of present polynomial value to last Taylor's series iteration polyno- 
mial value 

value of normalized polynomial at XHIGH 
value of normalized polynomial at XLOW 
square root term of quadratic equation 

square root term of Taylor’s series approximation for real root location 
polynomial ratio used to estimate location of next root 
temporarily real part of polynomial coefficient at a complex root 
radial component of a spiral in the complex plane 

ratio of a polynomial derivative to the next lower order polynomial deriva- 
tive 
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RATIOI 

RATIOR 

RRATIO(I) 

SCALE 

SLOPE 

SRATIO 

SXI 

SXIYI 

SXI2 

SYI 

THETA 

THETAO 

TOLI 

TOLIQ 

TOLR 

TOLRQ 

X 

XCON 

XCPLX 

XDT 

XH 

XHIGH 

XI 
XL 

XLAST 


imaginary part of polynomial ratio used to estimate location of the next root 
real part of polynomial ratio used to estimate location of next root 
real part of complex number derivative ratio 

hexadecimal order of magnitude factor used to scale polynomial roots 

slope of least squares line fit of the hexadecimal logarithms of the polyno- 
mial coefficients 

parameter used in setting step size when between or among roots of a 
group, V p OLY(ND + 1) * POLY(ND - 1) 

sum of XN for determination of SLOPE 

sum of the product of XN times YI 

SXI squared 

sum of YI for determination of SLOPE 
local angular coordinate of a spiral in the complex plane 
reference angle of a spiral in the complex plane 
complex root resolution distance criterion 

probable absolute error of a computed normalized polynomial value in the 
complex plane 

real root resolution distance criterion 

probable absolute error of a computed normalized polynomial value along 
the x-axis 

independent polynomial variable or real component of independent complex 
variable 

one of the constants used to determine DELTX and DELTY for a spiral 
x-coordinate at which the search for complex roots is begun 

real part of the denominator of the complex Taylor's series approximation 
for a new z 

temporary new value of x 

closest x definitely known to be greater than the value of the root sought 
initial x value in the search for a new root 

reference x when adjusting to a new x value between XHIGH and XLOW 
x value of a previous iteration 
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XLOW 

XMULT 

XN 

XNT 

XOLD 

XP(I) 

XPN 

XPOLE 

XPOLY(I) 

XQ 

XR(I) 

XROOT(I) 

XX 

Y 

YCON 

YCPLX 

YDT 

YI 

YLAST 

YNT 

YOLD 

YP(I) 

YPOLE 

YPOLY(I) 


closest x definitely known to be less than the value of the root sought 
number of roots at a root point 

number of a polynomial coefficient, which is the independent variable in 
the determination of SLOPE 

real part of the numerator of the complex Taylor’s series approximation 
for a new z 

x value of a previous iteration during resolution of a root 

real part of normalized complex polynomial or polynomial derivative 

temporary storage of x or XP(I) 

x-coordinate of spiral pole in the complex plane 

real part of complex polynomial or polynomial derivative 

real part of square root term in the complex number Taylor’s series solu- 
tion for a new z 

polynomial remainder or real part of polynomial remainder after a root is 
divided out 

x-coordinate of a root 
temporary storage of x 

imaginary component of independent complex variable 

one of the constants used to determine DELTX and DELTY for a spiral 

y-coordinate at which the search for complex roots is begun 

imaginary part of denominator of the complex number Taylor's series ap- 
proximation for a new z 

initial y value in the search for a new root, also hexadecimal logarithm 
of a polynomial coefficient 

y value of previous iteration on the spiral 

imaginary part of the numerator of the complex Taylor’s series approxi- 
mation for a new z 

y value of a previous iteration during resolution of a root 

imaginary part of normalized complex polynomial or polynomial derivative 

y-coordinate of spiral pole in the complex plane 

imaginary part of complex polynomial or polynomial derivative 
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on o o n o 


YQ 

YR(I) 

YROOT(I) 

YY 

ZCON 


imaginary part of square root term in the Taylor’s series approximation 
for a new z 

imaginary part of the polynomial remainder after a root is divided out 
y-coordinate of a root 
temporary storage of y 

one of the constants used to set up the spiral increment 


Program Listing 


$ I BF TC RCCTS 

DOUBLE PRECISION A, ARATIO, B, CA, CRATIO, D, D2, FK, P, POLY, Q, 

1 RA, RATIO, RRATIO, SCALE, X, XHIGH, XLAST , XLOW, XP, XPN, XPOLY, 

2 XR, XRCCT, Y, YP, YPOLY , YR, YROOT 

CCMMCN /FOLYN/ ARATIO, D, DELTX, DX, D2, FK, I, INK, IREAL, 

1 ITRIG, IW, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N, NC, ND, 

2 NCEG, NCRV , NEXP, NEXPS, NK, NM, NUM, OARAT, Q, SCALE, SRATIO, 

3 TCL I , TCLIQ, TOLR, TOLRQ, X, XI, XLAST, XMULT , XOLD, XPN, Y, 

A A ( ICO ) , B(100), CA(100I, CRATIO(IOO), P(IOO), POLY(IOO), RA(100), 
5 RATIC(ICO), RRATIO(IOO), XP(IOO), XPOLY(IOO), XR(100), XROOT(IOO) 
6 , Y P ( 100 ) , YPOLY ( 100 ) , YR(IOO), YROOT(IOO) 

ITEST = 1 

IF (ITEST. EO. 2) CALL DUBIO 

A ( I ) ARE COEFFICIENTS FOR THE POLYNOMIAL, POLY(X) = A(N)*X**N 

* A(N-1 )*X*=MN-1 ) + + A(2)*X**2 + A(1)*X + A(0). THE COEFFI- 

CIENTS MIST BE INPUT IN ORDER. START THEM FROM THE HIGH DEGREE 

ENC CF TFE POLYNOMIAL. 

IR = 5 

rv. = 6 

5 READ C I R , 1000 ) N, (A(I), 1*1, N) 

IF (N.LE.I) WRITE (IW,.1010) 

EF (N.GT.100) WRITE (TW,1020) 

WRITE (IV, 1030) (A(I),I*1,N) 

NUM = N -1 
J = C 
NC = NUM 

IF (NC.GT.10) NC = 9 + NUM/10 
TCLRC = 5.0*2.0**NC*1.0E-16 
TCL R =i SCRT(TOLRC) 

TCL I C = 2.0**NC*1.0E-15 
TC L I SCRT(TOLIG) 

NC * 1 

IF THE POLYNOMIAL BEGINS WITH ZERO COEFFICIENTS LOWER THE 

DEGREE APPROPRIATELY. 
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DC 6 1=1, N 

IF (MI ) .NE.O.O) GO TO 7 
NC 3 NC + 1 

6 CONTINUE 
GC TC 5 

7 IF (NC.EC.l) GO TO 9 
DC 8 I = N C , N 

NM - I - NC + 1 
E A(NK) = MI) 

N=N-f^C + l 
NUN - N - 1 
9 XI = 0.0 
NC EG = N- 1 

C TAKE CUT ANY ZERO ROOTS FIRST. 

FK = AB S ( A ( 1) ) 

DC 1C I =2 , NDEG 

IF <AeS<MI)).LT.FK*0.01*TOLRQ) GO TO 10 
FK = ABS< MI ) ) 

1C CONTINUE 

IF (ABS(MN)/FK) . GT . 0 . 0 1 *TOLRQ ) GO TO 12 
X = C.O 
J 3 J + 1 
XRCCT(J) = 0.0 
YRCCTfJ) 3 0.0 
XR ( J ) = C.O 
YR ( J ) = C.O 
N = N - 1 
GC TC 9 
12 IREAL = 1 

SCALE POLYNOMIAL SO ROOTS ARE APPROXIMATELY OF ORDER ONE FOR 

RCCT RESOLUTION COMPUTATION PURPOSES. A LEAST SQUARES FIT OF THE 

LCGS CF THE POLYNOMIAL COEFFICIENTS TO A LINE IS USED. 


12 SXI 3 O.C 
S X I 2 3 0.0 
SYI = O.C 
SXIYI = C.O 
NC 3 C 

FK 3 A B S ( A ( 1 ) ) 

DC 1 A I=1,N 

IF ( ABS ( A f I) ).LT .FK*0 .01*T0LRQ ) GO TO 14 
NC NC 4 1 
FK 3 A 8S < A ( I ) ) 

XN 3 FLOAT! I ) 

Y I 3 0.36067376* A LOG ( FK ) 

SXI ^ SXI + XN 
SXI2 * S> 1 2 + XN**2 
SYI 3 SYI + Y I 
SXIYI 3 SXIYI + XN*Y I 
14 CONTINUE 

SLCPE 3 (SYI/ FLO AT ( NC ) - SX I Y I / SX I) / ( SX I /FLO AT ( NC ) - SXI2/SXI) 
FK = 0. 36067376* ALOG(ARS(A( 1 ) ) ) 

NC 3 FK 
NK 3 C 
INK = 1 

EF (SLOPE .GT. 0.0 ) GO TO 15 
SLCPE = -SLOPE 
INK = -1 
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15 DC 16 1=1,10 

IF (SLOPE.LT. 0.5 ) GO TO 17 
NK = NK 4 INK 

16 SLCPE = SLCPE - 1.0 

17 NEXP = NK 
NEXPS = NEXP 

SCALE = 16.0D0**NEXP 
IF (NK.EC.O.AND.NC.EG.O) GO TO 19 
DC 18 1=1, N 
NM = NK* ( 1-1 ) + NC 
16 A ( I ) = A ( I )/16.000**NM 
IS IF (N.LE.3) GO TC 570 

IF (IREAL. EC. 0) GO TO 280 

FACT = 2.0 

ICIP = 1 

XCPLX = C.O 

YCPLX = 1.0 

INITIALISE PARAMETERS. X IS THE INDEPENDENT VARIABLE. J IS 

THE NUMBER OF ROOTS FOUND. FACT IS AN X STEP SIZE FACTOR. IDIR 

INDICATES THE DIRECTION TO MOVE ALONG X. IREAL IS A TRIGGER TO 

IN'CICATE THE SEARCH FOR REAL OR IMAGINARY ROOTS. ONLY REAL 

- — ALGEBRA IS USED WHEN IREAL IS A POSITIVE INTEGER. H IS A DELTA X 

INCREMENT. ITRIG IS A ROUTING DEVICE. IT BASICALLY GIVES AN 

INCICATICN OF HOW THE POLYNOMIAL APPROACHES THE X AXIS. JTRIG IS 

A ROUTING CEVICE WHEN HIGHER DERIVATIVES ARE NEEDEO. 

2C H = C.O 
DIR = IDIR 
IMSURE = 0 
ITRIG = C 
IXM = 0 
JTRIG = C 
KTRIG = C 
LTRIG = C 
X = XI 

ND IS THE NUMBER OF DERIVATIVES CALCULATED. 

NDRV IS A PARTICULAR DERIVATIVE. 

25 NC = 2 

THIS PROGRAM USES A QUOTIENT OF SUCCESSIVE DERIVATIVES METHOD 

— - IN BCTH THE REAL AND COMPLEX REGIMES TO PINPOINT THE VALUES AND 

THE MULTIPLICITY OF THE ROOTS. SINCE THE REAL ROOTS ARE EASIER TO 

WCRK WITH, THEY ARE FOUND FIRST. IF THERE ARE REAL ROOTS THE 

DEGREE OF THE POLYNOMIAL IS DECREASED BEFORE COMPLEX NUMBERS ARE 

USED. 

28 ICFECK = 0 
3C DC 32 1 = 1, N 
32 BID = A ( I ) 

NC EG = N-l 
M = NCEG 
NCR V = 0 
35 CALL RPCLY 

IF ( ICHECK.EQ.O.OR.AESl P ( 1 ) ) .GT .0. 5*T0LRQ ) GO TO AO 
IF ( ICHECK.GT.O) GO TO 38 
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36 X = tXLOW + X )/2 «0 
GC TC 28 

35 X = (XHIGH + X ) / 2 • 0 
GC TC 28 
4 C ICHECK = 0 

EF (ITRIG,EQ.O.OR.ABS(P( 1) ) . LE . 0 • 5+TOLRO ) GO TO 48 

IF ( IVSUPE.EQ.O) GO TO 44 

IF (P(l )/PXHIGH*LT.O.O) GO TO 42 

XHIGH = X 

PXHIGH = P(l) 

GC TC 48 
42 XLCU = X 

PXLCV = F 1 1 ) 

GC TC 48 

44 EF < PCLY<1)/0P0LY1.GT,0.0> GC TO 48 
IFSURE — 1 

IF ICELTX .LT.0.0 ) GO TO 46 

XHIGH = X 

PXHIGH = P ( 1 1 

XLCR =; XL AST 

PXLCW = C PCLY1/B ( N ) 

GC TC 48 

46 XHIGH = XLAST 

PXHIGH = OPOLY1/ BIN) 

XLCW - X 
PXLCW = F(l) 

48 IF (AES(F(l)).LT.TOLR.OR.JTRIG.EQ.2) GO TO 160 
IF (IREAl.EO.O) GO TO 50 

IF ( ITRIC.NE.O.OR. ABS(P< 2) l.GT.TOLR) GO TO 50 
H = 0.1+CIR 

IF (ABS(FOLY( D) .GE . ABS ( POLY ( 3 ) ) ) GO TO 49 

IF ( POLY ( 1 ) / POLY (3) .LT.0.0) H = -D IR/ 1 POLY (1) *POLY ( 3 ) ) 

45 DELTX = H 
DIR = l.C 
ITRIG = -1 
GC TC 102 

C REAL ROOTS ARE THE X VALUES WHICH MAKE P(l) = 0.0 START AT 

C X = C.O. THE NEXT X IS DETERMINED WITH A 3 TERM TAYLORS SERIES. 

5C Q1 » (P0LY(2)/P0LY(3) )**2 - 2.0*P0LY (1) /POLY< 3) 

C IF Cl IS LESS THAN ZERO, P(2) MAY CHANGE SIGN BEFORE P(l) 

C CRCSSBS THE AXIS. 

IF (Cl. LT.0.0. OR. ITRIG. EQ.l) GO TO 70 
IF (JTRIG.EC.O) GO TO 57 
IF (PCLY(1)/0P0LYI. LT.0.0) GO TO 56 
IF ( ITRIG. LT.O) GO TO 52 
IF ( POLY ( 2 ) /POLY ( 1 ) *H ) 57,57,66 
52 IF ( POLY ( 2 ) /OPOL Y2 .GE .0 .0 ) GO TO 57 

56 DIR * l.C 

57 ITRIG = -1.0 

58 IF (CABS(P(2) l.LT.TOLRQ) P0LY(2) = T0LRQ*CA<2) 

IF ( ABS ( FOLY (3)/ POLY ( 2) ) .GT .0.01 ) GO TO 60 

55 H = — POLY ( 1 ) /POL Y(2)*DIR 
IF (X) 64,65,64 
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6C IF (ABS(F (3) ) .LT.TOLRQ) GO TO 62 
IF (ABS(F(2)) .LT .TOLRQ) GO TO 61 

H * PCL Y ( 2 ) /POLY (3)*(SQRT(Ql)/ABS(POLY(2)/POLY(3)) - 1.0) 

GC TC 63 

61 H = CIR+SQRT (Ql) 

GG TC 63 

62 H = - PCL Y ( l ) /POLY ( 2 ) 

63 IF (ABS(P) . LT. 100. 0*T0LR. AND. ABS( POLY t 2)/POLY(l) ) .GT. . 01 /TOLR ) 
X GC TC 55 

IF (ABS(P ).LT.1.00) GO TO 65 
IF (ABS(X).GT.O.l) GO TO 64 
H » H/ABS ( H ) *0. 1 
GC TC 65 

64 IF (H/X.GT.0.5) H = 0.5*X 
IF (H/X.LT.-1.5) H * -1 . 5*X 

65 DELTX = P *C I R 
GC TC 68 

66 IF (ITRIG. GE. 15) H = 1.5*H 
DELTX = P 

ITRIG = ITRTG + 1 
68 JTRIG = -1 
GC TC 102 

7C IF (IREAL. EC. 0) GO TO 275 

C STATEMENTS 70 THROUGH 77 INVESTIGATE Ql LESS THE 0. 

72 IF ( ITRIG. GT. 15) GO TO 110 

IF (ITRIG. GT. 10. AND. ABStP(l) ) . GT . 1 .O/TOLRQ ) GO TO 110 
IF (JTRIG.GT.O.OR.ITRIG.GE.2) GO TO 74 
HI = -PCL Y ( 2 ) /POLY ( 3 )*DI R 

IF (ITRIG. EQ.O. AND. ABS(H1).GT. 1.0) HI * H1/ABS(H1) 

H = PI 

IF (X.EQ.XI) GO TC 73 

IF (ABS(F1/X).GT .1.0) HI = H1*ABS(X/H1) 

H * FI 

IF ((X -X I)*DIR/'H.LT.O.O) H = -H 

73 TF ( ITRIG. LE.O) ITRIG = 1 
JTRIG = 1 

PCLY1R = 1.0 
GC TC ICC 

74 PCLY1R = PCLY(l) /0P0LY1 
TF (JTRIG. LT.O) HI = H 

IF (PCLY1R.LE.0.0): GO To 75 

IF (PCLY(2)/0P0LY2.LT. 0.0. AND.ITRIG.LT. 2) GO TO 76 

C IF F0LY1R IS GREATER THAN 0.5. NO ROOT IS EXPECTED IN THIS 

C REGION CF Cl LESS THAN 0. 

IF (P0LY1R.GT.0.95) GO TO 80 
Hl= F1*FACT 
GC TC 78 

75 ITRIG = -1 

IF (FACT. LT. 1.0) GO TO 77 
FACT s 0.5 
H = HI 
H = -F*F A CT 
GC TC 90 
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THE POINT OF 


THE SIGN OF THE SLOPE OF POLY(X) HAS CHANGED. 

INTEREST HAS BYPASSED. HALF STEP BACK. 

76 IF (ABS(H).LT.TCLR) GO TO 80 

77 FACT* 0.5 
HI* -H1*F ACT 

76 XCPLX = > 

YCPLX = SQRT ( ABS ( POLY f 1 ) ) ) 

GC TC IOC 

IF STATEMENT 80 IS REACHED, POLY(X) DID NOT CROSS THE X AXIS. 

CONTINUE IN THE SAME DIRECTION THAT WAS STARTED EVEN THOUGH P(X) 

IS MOVING AWAY FROM THE X AXIS. THE STEP SIZE IS DOUBLED UNTIL 

THE SECOND DERIVATIVE CHANGES SIGN. 

8C IF (IREAL.EQ.O) GO TO 275 
IF (ITRIG.LT. 2) H = O.l+H 
FACT ^ 2.0 
ITRIG * ITRIG + 1 
H = h*2 • C 

SC IF (H. EC. 0.0) H = 0.1+DIR 
DELTX = F 
GG TC 102 
IOC DELTX * FI 
102 0PCLY2 = POL Y ( 2 ) 

IF (GP0LY2.EQ.0.01 0P0LY2 = 1.0E-17 
10^ IF (ABS(X).GT.l.OI GG TO 106 

IF (ABS(CELTX).LT.TOLRQ) GO TO 160 
GC TC 108 

106 IF (4BS(CELTX/X) .LT.TOLRO) GO TO 160 
108 GPCLY1 * PCLY(l) 

XL AST = X 

X * X + CELTX 

IT I I MSUP E. EO *0 ) GO TO 30 

IF IX.LE.XHIGH.AND.X.GE.XLOW) GO TO 30 

X = tXHIGH + XLOW ) /2 • 0 

HI (XHIGH * XL0WJ/2.0 

IF (XHIGF.EO.XLAST) HI = -HI 

DELTX = FI 

GC TC 30 

STATEMENT 110 IS REACHED WHEN THERE ARE NO REAL ROOTS FOUND 

IN THE DIRECTION STARTED. GO BACK TO XI AND WORK THE OTHER 

DIRECTION IF IT HAS NOT BEEN INVESTIGATED. 

1IC IF (INSURE. EG. 0) GO TO 115 
DELTX * (XHIGH - XL0W)/2.0 
IF (X. EC. XHIGH) DELTX = -DELTX 
GC TC 102 

115 IF (ICIR) 130*130* 120 
12C IDIR = -1.0 
X = XI 
GC TC 20 

STATEMENT 130 IS REACHED IF NO REAL ROOTS WERE FOUND IN 

EITHER SIDE OF X ^ XI. SEARCH FOR COMPLEX ROOTS. 
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13C IREAL = C 

X = XCPLX 
Y = YCPLX 
GC TC 28C 

14C IF (ABS(F(1> ).LT.0.5*T0LRQ) GO TO 142 
OPCLY1 = PCLY(l) 

XLAST = X 

IF ( ITRIG . EC. 0) ITRIG = -1 
142 IT (IMSUFE.EQ.O) GO TO 150 
XH = X + DELTX 

IF (XH.LT. XHIGH. AND. XH.GT. XLOW) GO TO 150 
IF (LTRIG.EC.l) GO TO 193 
IF (KTRIG.EQ.2) GO TO 1244 
IF (LTRIG.LT. 0) LTRIG = LTRIG - 10 

IF (ABS(FXHIGH).LT.TOLRQ.ANO.ABS(PXLOW).LT. TOLRQ ) GO TO 155 
XL = 0.0 

IF (CELTX.GT.0.0 ) GO TO 145 
14 2 TF (AeS(FXLCW).LT. TOLRQ) GO TO 145 
IF (X.NE.XLOW) GO TO 144 
XL = -0.5*DELTX 

IF (XHIGH .GE.X-OELTX) GO TO 144 
XL = SORT (A8S IPX LOW ) ) 

XI * (XHIGH - XLOW ) *XL/ ( XL + SQRT( ABS( PXHIGH) ) ) 

144 DELTX = XL + (XLOW - X - XL)/t0.7 +0. 5*AL0G ( ABS ( PXL0W*2 . O/TOLRQ) ) ) 
I CHECK = -1 

GO TC 146 

145 IF (ABS(FXHIGH).LT. TOLRQ) GO TO 143 
IF (X.NE. XHIGH) GO TO 146 

XL => — 0 .5 *DELTX 

IF (XLOW . LE .X— DELTX ) GO TO 146 
XL * SORT (ABS(PXHIGH) ); 

XL = (XLCW - XHI GH) *XL/ ( XL + SORT ( ABS ( PXLOW )) ) 

146 DELTX = XL + (XHIGH -X -XL)/(0.7 +0. 5* ALOG ( AB S ( PXHIGH*2 . O/TOLRQ )) ) 
ICHECK * 1 


148 

CF 

(LTRIG. LT.O) 

GO TO 

150 


KTR 

IG 

= C 




IX* 

= 

0 




NC 

* 

2 



14 5 

LTR 

IG 

= C 



15C 

X = 

X 

+ CELTX 




HI 

a 

CELT X 




GC 

TC 

30 



152 

V = 

F 

- 1 




DC 

15 

4 1*1, M 




FK 

* 

N - NCRV - 

I + 1 


154 

B< I 

) 

= B ( I ) *FK 




GC 

TC 

35 



155 

IF 

(LTRIC.LT. 0) 

GO TC 

201 

158 

X = 

(XHIGH + XLOWJ/2.0 



GC TC 251 

C STATEMENTS 160 THROUGH 250 ARE THE RATIO OF DERIVATIVES 

C METHCC FCR REAL ROOTS. 

16C IF (ITRIG. GE. 2) GO TO 80 

IF ( ABS (F ( 1 ) ).GT. TOLRQ. AND. ABS(POLY( 2)/P0LY(l) ) .LT. 1.0) GO TO 165 
IF (NC. EC. NDEG. AND. LTRIG. EQ.O) GO TO 170 
IF (ABS(F lND-1)) .LT. TOLRQ) GO TO 162 
LIMIT = C 
GC TC 18C 
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180 


162 UNIT = 1 

IF (JTRIG.LT. 2) GO TO 164 
IF (ABS(F(ND) r.GT.TOLRQ) GO TO 
DIN IT = 4 
GC TG 251 

164 IF (ABS(F(2) ).GT.10.0*T0LR) GO TO 180 

ND =? ND 4 1 

GO TC 152 

165 QI = ( POL Y ( 2 ) /POLY ( 3 ) V**2 - 2 .O+POLY ( 1 ) /POL Y( 3 ) 

166 JTRIG = 1 
KTRIG = C 
LTRIG = C 
NC = 2 
xcplx = y 

YCPLX = C.001 
IF (Cl) 167,168, 168 

167 ITRIG = I 
GO TC 80 

166 IREAL = 1 
GC TC 58 

17C IF ( ABS ( F ( ND ) ).L E .0. 5*T0LR0 ) GO TO 252 
RATICCND) = POLY ( ND+ It/ POLY ( ND ) 

LIMIT = C 

IF ( AB$(PATIO(ND) ) . LT .0 . 5/TOLR ) GO TO 175 
TF (ABS(F (ND-1>> .LT.TOLRQ) GO TO 176 
RATIC(ND-I) = P0LY(ND)/P0LY(ND-1 ) 

IT (ABS(PATIO( ND ) /R AT 10 ( ND-1 ) ).GE.1.0) GO TO 180 
NC * NO 4 1 
GC TC 245 

175 IF (ABS(F<ND-1)) .GT.0.5*TOLR0 ) GO TO 180 

176 DELTX = TOLR 
GC TC 14C 

1 8 C IF (JTRIG. GE. 2) GO TC 181 
JTRIG = 3 
CALL SCALER 

IF (NK.EC.O ) GO TO 181 

X = X/16.0C0**NK 

ITRIG = C 

JTRIG = 2 

INSURE = 0 

GC TC 30 

181 IF ( ABS ( F ( ND-1 ) ) .LT .0 *5*T0LRQ ) POLY(ND-l) = 0. 5*T0LRQ*CA ( ND-1 ) 

182 RATIC(NO-l) = POLY ( ND >/ POLY ( ND- 1 ) 

185 RATIC(ND) = POLY ( ND+lt / POLY ( ND ) 

DRATIC = RATIO(ND)/RATIO(N D— 1 ) 

ARATIC = ABS(DRATIO) 

NEAP A ROOT WITH MULTIPLICITY EQUAL TO OR GREATER THAN ND -1, 

THE ABSOLUTE VALUE OF ARATIO SHOULD BE LESS THAN ONE. IF ARATIO 

IS GREATER THAN ONE, EITHER A ROUNDOFF ERROR IS ENCOUNTERED OR WE 

ARE NGT RELATIVELY CLOSE TO A ROOT. IN EITHER CASE USE RATI 0 ( ND ) 

TO MOVE >. HOWEVER, LIMIT THE MAGNITUDE OF RATIO(ND) SO THAT THE 

MULTIPLICITY CHECK CAN BE MADE WITHIN COMPUTER ACCURACY. 

IF (KTRIC.EQ.il GO TC 240 
IF (LTRIG) 200,186,190 

186 IF (ARATIO.LT. 1.0) GO TO 210 

186 SR AT I C = SCRT(ABS(RATI0(ND)*RATI0(ND-1) ) ) 

IF (SRATIO.LT.l.O) SRATIO = 1.0 
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19C IF (LTRIG.EQ.l) GO TO 192 
IF (LTRIG. GT.l) GO TO 1195 
XX * X 
IXN * 1 
JTRIG = 2 
LTRIG = 1 
OCR AT = CRATIO 
DELTX = C.l/SRATIO*DIR 
CTX * CELT X 
GO TC 1 AC 

192 IF ( CP AT IO/ODRAT .LT.O.O) GO TO 194 
IF ( ARATIO.LT. ABS(ODRAT) ) GO TO 196 

193 X = X - CELTX 
DELTX = -DELTX 
DX * DELTX 

GC TC 196 

194 DELTX = -DELTX*ARATIC/( ARATIC + ABSIODRAT ) ) 

OCR AT = CRATIO 

GO TC 14C 

195 IF (CRATIO/ODRAT. LT.O.O) GO TO 194 
IF (LTRIG. LT.O) GO TC 196 

1195 IF ( ARATIO.GT.ABS (ODRAT ) ) GO TO 199 

196 IT (ARATI0.LE.0.45) GO TO 230 

TF ( LTRIG.EC.8.AND. INSURE. EQ.O) GO TO 272 

197 DELTX * 2.0*DELTX 
OCR AT * CRATIO 

196 IF (LTRIG. GE.O) LTRIG = LTRIG + 1 
GC TC 14C 

199 IF ( INSUPE.EQ.O) GO TO 276 
LTRIG = -1 
IXN * 0 
GC TC 197 

2CC IF (LTRIG .GE. -10 ) GO TO 195 

201 IF (LTRIG. LT. -11 ) GO TO 202 
LTRIG = -2 

X = XX 
DELTX = -DX 
GC TC 14C 

202 IF (LTRIG.LT. -12) GO TO 204 
202 IF (AC. EC. 2) GO TO 158 

ND » NO - 1 

DX = (XH IGF - XL CW ) /10.0 
DELTX = CX 
X * XLOVf + DX 
LTRIG = -13 
GC TC 30 

204 IF (ARATIO.LT. 0.45) GO TO 230 

205 IF (LTRIG. LT. -13 ) GO TO 206 
GC TC 207 

206 IF (LIMIT. GT.O) GO TO 207 

IF (CRATIO/ODRAT. LT.O.O. OR. LTRIG. EQ. -15) GO TO 209 

207 X = X + CX 
LTRIG = -14 

IF (X.GE.XFIGH) GO TO 203 
206 CDRAT = CRATIO 
GC TC 30 
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2CS OELTX = C ELTX*DR ATIO/(ODRAT - ORATIOI 
LTRIG = -15 
CCRAT = CRATIO 
GC TC 15C 

21C IF (JTRIG. EQ. 4) GO TO 212 

IF (ARATIO - 0.45) 230,230,215 
212 JTRIG = 2 
NC =* ND - 1 
GC TC 30 

215 IF (LIMIT. EC. 0. OR. JTRIG. EQ. 3) GO TO 220 
GC TC 186 

216 JTRIG = A 

22C IF (AC.EC.NDEG) GO TO 225 
NC * NO ♦ 1 
GC TC 152 

225 LTRIG = C 

DELTX = -l.O/RATIOIND) 

IF (JTRIG. EQ. 4) CELTX - -1 .O/RAT 10 ( ND-1 ) 

JTRIG = 2 
GC TC 14C 

226 IF (LIMIT. EQ.l) GO TO 248 
IF (LTRIG. LT.O) GO TO 229 

22 6 NC 58 NO - 1 
GC TC 30 
225 DELTX = CX 
GC TC 207 

RATIC(ND-l) IN STATEMENT 230 IS THE LAST DERIVATIVE RATIO 

THAT SHOLLC APPROACH INFINITY AS X APPROACHES A ROOT VALUE. 

THE PROXIMITY OF X TO THE ROOT IS 1 .O/RAT 10 ( NO- 1 ) . 

THE MULTIPLICITY OF ROOTS AND THE ESTIMATE OF THE NEXT ROOTS 

ARE MADE WHEN THE POLYNOMIAL DERIVATIVE RATIO LIES BETWEEN 0.5T0LR 

ANC 5.0TCLR. THIS MEANS THAT ROOTS CLOSER THAN TOLR TOGETHER 

CANNCT BE RESOLVED FROM ONE ANOTHER. THE ACTUAL X VALUE OF THE 

RCCT , HOWEVER* IS IMPROVED SOMEWHAT IN STATEMENT 250. 

23C IF (LTRIG. GT.O) LTRIG = 0 
IF (AD. EC. 2) JTRIG = 2 
IF (JTRIG. EC. 3) GO TO 218 

IF ( ABS(RATI0(ND-1 ) ) .LT.5.0/T0LR ) GO TO 235 

232 IF ( ABS(FATIOlND-l) ) . LT . 1 . O/TOLRQ ) GO TO 233 
DELTX = 1.0/ (T0LR*RATI0(ND-1 )) 

GC TC 15C 

233 DELTX = -0.2*T0LR 
GC TC 15C 

235 IF ( ABS(FATIG(ND-1) > .GT . 1 .O/TOLR ) GO TO 242 
IF ( LIMIT. GT.O) GO TO 238 

DELTX = -1.0/RATIO(ND-1) 

GC TC 14C 

236 IF (AeS(PATIO(ND-D).LT. 1000.0) GO TO 251 
XCLC = X 

CRATIC = RATIO(NC-l) 

K T R T C = 1 

DELTX = C.l/ORATIO 
CRATIC = SCRT( ABS(ORATIO) ) 

GC TC 15C 
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24C IF ( ARATIO. GE.O. 45) GO TO 241 

IT ( ABS ( RAT 10 (ND-1 ) ) .LT.ORAT 10 ) GO TO 241 
IT (LIMIT. EG. 0) GO TO 248 
DELTX = 2.0*DELTX 
GO TC 14C 

241 X = XOLD 
GC TC 251 

242 IF (NC.LE.2) GO TO 245 

IF (LIMIT. GT.O) GO TO 1245 

IF ( ABS ( P ( NO-2) ) .GE.0.5*TOLRQ) GO TO 243 

LIMIT = 1 

PCLY (ND-2 ) = 0.5 *TOLRQ*C A ( ND-2 ) 

243 RATIO (NO-2) = POLY ( ND-1 ) /POL Y( NO-2 ) 

ARATIC = ABS( RAT 10 (ND-1 ) /RAT I0( ND-2) ) 

IF (KTRIG. EC. 2) GO TO 244 

IF (ARATIO.LT. 0.53) GO TO 248 

IF (IXM.EQ.O) GO TO 226 

CARAT = ARATIO 

XPN = X 

KTRIG = 2 

GC TC 247 

244 IF (ARATIO.LT .0. 53) GO TO 248 

IF (ARATIO. GT.ABS(ODRAT) ) GO TO 1244 
IF (LIMIT. EC. 0) GO TO 228 
GO TC 246 

1244 X = XPN 
KTRIG = C 
IXM = 0 
DELTX = C.C 
GC TC 14? 

C THE MULTIPLICITY OF THE ROOT AND ITS VALUE WITHIN TOLR ARE 

C NOW AT HANC. ESTIMATE THE X LOCATION OF THE NEXT ROOT. 

245 IT ( ABS(F0LY(3) ) .GT.ABS( POLY ( 1 ) ) ) GO TO 248 

1245 IF (IXM.EQ.O) GO TO 248 

246 I"XM = 0 

247 DELTX = 2.0*(XX - X) 

LTRIG = C 

GC TC 14C 
246 MULT =; NC - 1 
XMLLT = MULT 

IF (NCEG - MULT.LT.3) GO TO 250 
RA = 1.0/ ( XMULT + 1.0)*RATIO(ND) 

XI = X -1.0/RA 
IT (ABS(RA).LT.O.l) XI = X 
25C X = X - 1.0/RATI0(MULT) 

GC TC 255 

251 NO = 1 

252 MULT =: NC 
XI a X 

C NUMBER THE ROOTS AND DIVIDE THEM OUT OF THE POLYNOMIAL. 

255 DC 270 K=1,MULT 
J = J +1 

XRCCT(J) m X*SCAL E 
VRCCT(J) = 0.0 
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N = N— 1 

IF (N.EC.l) GO TO 265 
DC 260 1*2, N 

26C MI) = All) + X*AII-1) 

265 XR(J) * 1.0 + X* A I N ) / At N + l ) 

27C YRIJ) = C.O 
GC TC 19 
272 X = XX 

GC TC 276 

275 IREAL * -1 

276 X = X — CFLTX 

276 Y * SORT ( ABS ( POL Yl 1)1): 

28C CALL COMPLX 

TF (IREAL.LT. 2) GO TO 19 
IREAL = C 
GC TC 252 
57C CALL CEO 

10CC FORMAT ( 1 3/ ( 3D24 . 16 ) ) 

101C FORMAT (//25X*78HERR0R THE DEGREE OF THE POLYNOMIAL IS LESS T 

XHAN CNE. THERE ARE NO ROOTS. ) 

102C FORMAT (//37X*,55HERR0R THE PROGRAM IS ONLY DIMENSIONED FOR N = 

X ICC. ) 

103C FORMAT < 1H1 ft //3 IX , 70HT HE INPUT POLYNOMIAL COEFFICIENTS. THEY ARE 
X IN CRDEP IF READ IN ROWS. U (2X,5D24.16>) 

GO TC 5 
ENC 


SLERCUT IN E RPOLY 

THIS IS THE BASIC ROUTINE FOR FINDING THE VALLE OF A REAL 

NUMBER POLYNOMIAL AND ITS DERIVATIVES. P( ) ARE THE POLYNOMIAL 

ANC ITS DERIVATIVE VALUES NORMALIZED BY BIN). 

DOUBLE PRECISION A, ARATIO, B, CA, CRATIO, D, D2, FK, P, POLY, 0, 

X RA, RATIO, RRATIO, SCALE, X, XL AST , XP, XPN, XPOLY, XR, 

X XRCCT, Y , YP, Y POLY, YR, YROOT 

COMMON /FOLYN/ ARATIO, D, DELTX, DX, D2, FK, I, INK, IREAL, 

1 ITRIG, IW, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N*, NC , ND, 

2 NCEG , NCRV, NEXP, NEXPS, NK, NM, NUM, OARAT, Q, SCALE, SRATIO, 

3 TCLI, TOL IQ, TOLR , TOLRQ, X, XI, XL AST , XMULT , XOLD, XPN, Y, 

A A ( ICO ) , B(100), CA(IOO), CRATIO(IOO), P(100), POLY(IOO), RA(IOO), 

5 RATIO! ICO ) , RRATIO(IOO), XP(IOO), XPOLY(IOO), XR(IOO), XROOT(IOO) 

6 »YP ( ICO ) , YPOLY(IOO), YR(IOO), YROOT(IOO) 

35 PtNCRV+l ) * 1.0 
L = M 

TF (M.EC.O) GO TO 45 
NC = NEXFS - NEXP 
IF (L.EC.N-1) GO TO 37 
NM * NC 

36 IF (ABS(E(L+1I/CA(N0RV) J.GT. 0.001*8. 0**NM*T0LRQ) GO TO 37 
L * L - 1 

IF (L.EC.O) GO TO 37 
NM => NM ■* NC 
GC TC 36 
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37 I * 1 

36 IF (I.EC.L+1) GO TO 42 

PINCRV+l ) = P INDRV+1 ) *B I I ) 

II - I 
NP = NC 

35 IF (AeS(E(I+l)/B(Iin.GT.O. 001*8. 0**NM*T0LRQ) GO TO 40 
P (NCRV+1 ) = P (NDRV+1 ) *X 
IF (I.EC.L) GO TO 41 
1 = 1 + 1 
NP = NM + NC 
GO TC 39 

4C PINCRV+l) = P!NDRV+1)*X/B( 1+1) + 1.0 
1 = 1 + 1 
GO TC 38 

41 PINCRV+l) = P(NDRV+1)/B( I + 1 ) + 1.0 

42 IF (L.EC.M) GO TO 45 
PINCRV+l) = P(NDRV+1)*X**(M— L) 

45 POLY (NDRV+1 ) =i P ( NDRV + 1 ) *B I L +1 ) 

CAI NCRV + 1) = BIL + l) 

NCRV =! NC R V + 1 
IF (NCRV. GT. NO) RETURN 

48 P = P-1 

DC 49 1=1, P 

FK = N -NDRV -I +1 

49 BII) B(I) *FW 
GC TC 35 

ENC 


SUERCUT IN E COMPLX 

C THIS SUBROUTINE CONTAINS THE ANALYSIS PROCEDURES FOR COMPLEX 

C RCCTS. 

DCL6LB PRECISION A, ARATIO, B, CA, CRATIO, D, D2, FK, P, POLY, Q, 

X RA, RATIO, RRATIO, SCALE, X, XLAST , XP, XPN, XPOLY, XR, 

X XRCCT, Y, YP *. YPOLY, YR, YRCOT 

CCPPCN /FOLYN/ ARATIO, D, DELTX, OX, 02, FK, I, INK, IREAL, 

1 ITRIG, IW, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N, NC , ND, 

2 NCEG , NCRV, NEXP, NEXPS, NK , NM, NUM, OARAT, Q, SCALE, SRATIO, 

3 TCL I , TCLIO, TOLR , TOLRQ, X, XI, XLAST, XMULT , XOLD, XPN, Y, 

4 A ( ICO ) , BI100), CA(IOO), CRATIO(IOO), PI100), POLY(IOO), RA(IOO), 

5 RATIO! ICO) , RRATIO ( 100 ) , XP(IOO), XPOLY(IOO), XR(IOO), XR00TI100) 

6 »YP ( 100 ) , YPOLYI 100 ) , YR(IOO), YROOT(IOO) 

THE FIRST X ESTIMATE FOR THE FIRST COMPLEX ROOT IS SOME LOCAL 

MINIMUM CF THE ABSOLUTE VALUE OF POLY(X). LET THE 1ST Y ESTIMATE 

ECUAL THE SOUARE ROOT OF POLY(X). 

DTFETA = 3.1415927/6.0 
28C IF (IREAL. EC. -1) GO TO 528 
IF (Y.LT. 0.001) Y = 0.001 
IF (Y.GT.1C.0) Y = 10.0 
IF (AESIX l.LT.0.01) X = 0.01 
NCEG = N - 1 
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29C JTRIG = 1 
ITRIG = C 
IXM = 0 
NC * 2 

295 LTRIG = C 
30C DC 310 I * 1 »N 
31C B ( I ) = All) 

H- NCEG 
NCRV= 0 

32C CALL CPOLY 

351 IF (JTRIG. GE. 2) GO TO 360 

IF (ABS(XP(1) ).LT.TOLI.AND.ABS(YP( 1) ).LT.TOLI ) GO TO 360 

NEW X AND Y VALUES ARE DETERMINED WITH A 3 TERM COMPLEX 

- — TAYLCRS SERIES. 

XC => IXPCLY(2)/B(N) )*#2 - ( YPOLY ( 2 ) /B( N ) )**2 + 2.0*(YP(1)* 

X YPCLY ( 3 ) - XP(1 )*XP0LY(3) )/B(N) 

YC * 2.0* (XPOLY( 2)*YPOLY(2)/B(N) - ( XP ( 1 ) *YPOLY ( 3 J + YPC1)* 

X XPCLY < 3)) ) /B(N) 

PHI = ATA N2 ( YQ»X C ) 

IF (PHI. LT. 0.0) PHI = 6.2831853 + PHI 
PHI = PHI/2.0 

HYP - ( XC * *2 + Y Q**2 ) **0 .25 
XC = HYP*CCS( PHI ) 

YC = HYP*SIN(PHI ) 

DXT = X PC LY ( 2 ) /B ( N ) 

DYT = YPC LY ( 2 )/B ( N ) 

DZT = (XC - DXT)**2 + (YQ - DYT)**2 
DZT2 =5 ( > 0 + DXT ) **2 + (YQ + DYT)**2 
IF (CZT.LE.0ZT2) GO TO 352 
XC - -XC 
YC = -YC 
DZT = DZT 2 

352 XNT = XC - DXT 
YNT ^ YC - DYT 

XCT = X PC LY ( 3 ) /B ( N ) 

YDT = YPCLY(3)/B(N) 

D = XCT**2 + YDT **2 

HX s (XNT*XDT ♦ YNT*YDT)/D 

HY - (YNT *XDT - XNT *YDT )/D 

IF (ABS (HX/X) .LT .0. 1*T0L I.AND.ABS(HY/Y).LT.O. 1*T0L I ) GO TO 360 
TF (AES(HX).LT.l.O) GO TO 357 
IF (ABSO).GT.O.l) GO TO 356 
HX = HX / A BS ( HX ) *0 • 1 
GC TC 357 

356 IF (AES(FX/X).GT .0.5) HX = 0 .5*HX*ABS ( X/HX ) 

357 IF (ABS(HY).LT.l.O) GO TO 359 
IF (ABS(Y).GT.O.l) GO TO 358 
HY = HY/ABS(HY)*0.1 

GC TC 359 

356 IF (ABS(FY/Y) .GT.0.5) HY = 0 .5*HY*ABS( Y/HY ) 

359 X = X + HX 
Y= Y+HY 
GC TC 30C 

STATEMENTS 360 THROUGH 535 COVER THE RATIO OF DERIVATIVES 

METHCC FCR COMPLEX ROOTS. 
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36C IF ( ASS ( XP (ND-1 ) ) .LT. 0. 5*T0L IQ. AND. ABS ( YP ( ND-1 ) ) .LT.O. 5*T0LIQ) 

X GC TC 370 
LIMIT = C 
GC TC 38C 

c WHEN LIMIT = It POLYIND-1) HAS ROUNDOFF LIMITATIONS. 

37C LIMIT = 1 

IF (JTRIG. EQ.l) GO TO 375 

C WHEN LIMIT = 4, POLY(ND) HAS GREATER ROUNDOFF LIMITATIONS. 

rF (ABS(XP(ND)).GT.0.5*T0LIQ.OR.ABS(YP(ND) ) . GT. . 5*T0L IQ ) GO TO 380 
LIMIT = 4 
GC TC 512 

375 IF (AC.GT.NDEG/2 ) GO TO 380 
NC = NO 4 1 
GC TC 482 

38C IF (JTRIG. EC. 2) GO TO 382 
CALL SCALER 

IF (NK.EC.O) GO TO 382 
X = X/16.0C0**NK 
Y = Y/16.0C0*#NK 
JTRIG = 2 
GC TC 30C 

382 IF ( ABS(XP(ND-1) ) . LT .0. 5 *TOL IQ ) XP(ND-l) * 0.5*T0LIQ 
IF (ABS (YP(ND-l) ).LT.0.5*T0LIQ) YP(ND-l) = 0.5*T0LIQ 
D = {XP(f^D-l ) )**2 + ( YP ( NO- 1 ) ) **2 
PCLY(ND-l) * DSQPT(D)*DABS(CA(ND-1) ) 

IF ( ABS(XP(ND)).LT.0.5*T0LIQ) XP(ND) = 0.5*TOLIQ 
IF (A8S(TP(N0n.LT.0.5*T0LIQ) YP(ND) = 0.5*TOLI0 
D1 * (XP(NC))**2 4- (YP(ND) )**2 
PCLY(ND) = DSQRT (D1)*DABS(CA(ND) ) 

RAT IC (ND- 1 ) = POLY(NO)/POLY(ND-1) 

IF (ND.LE.2.0R.L IMIT.GT.O) GC TO 385 

IF ( ABS ( X P ( ND-2 ) ) .LT.O. 5*T0L IQ ) XPIND-2) = 0.5*T0LIQ 

IF ( ABS ( Y P (ND-2 ) ) .LT.O. 5*T0L IQ ) YP(ND-2) = 0.5*T0LIQ 

PCLY ( ND-2 ) * DSQRT ( ( XP ( ND-2 ) )**2 4- ( YP ( ND-2 )) **2 ) *DABS ( CA ( ND-2 I) 

RAT IC ( ND- 2 ) = P0LY(ND-l)/P0LY(ND-2) 

IF (LTRIG.NE.O) GO TO 385 
ARATIC = RATI0(NC-l)/RATI0(NC-2) 

IF ( ARATIO.LT. 0.45) GO TO 440 

IF ( ABS (XP (ND-2) ) . LT . 0. 6*T0L IQ . AND . ABS ( YP ( ND-2 ) ) .LT .0. 6*T0L I Q ) 

X GC TC 38 5 

IF < ARATIG.GT.1.0) GO TO 455 

385 PCLY (ND+1 ) * DSQRT ( (XP0LY(ND+1 ) )**2 + ( YPOL Y ( ND+ 1 ) ) **2 ) 

RATICtND) = POLY (ND+1 )/POLY(ND) 

ARATIC = RATI0(ND)/RATI0(ND-1) 

RRATIC(Nt-l) = ( XP ( NC ) *XP( ND-1 )+YP(ND)*YP( ND-1) ) /D*C A ( ND ) /CA ( ND-1 ) 
CRATIO(NC-l) = ( YP(ND) *XP( ND-1) -XP(ND)*YP (ND-1) ) /D*C A ( ND ) /C A ( ND-1 ) 
D3 = (RR f T I C ( ND- 1 ) )**2 + ( CR ATI0(ND-1 ) )**2 

RRATIC (NC ) = (XP (ND+1 )*XP( ND )+YP( ND+1) *YP(ND) ) /D1*C A ( ND+ 1 ) /C A ( ND ) 
CRATIC(NC) = (YP (ND+1)*XP(ND) -XP ( ND+1 ) *YP ( ND ) ) /D1*C A ( ND+ 1 ) /C A ( ND ) 
02= (RRA7 IC(ND ) ) **2 +(CRATIO(ND) )**2 
JTRIG = 2 

IF (LTRIG) 389,387,395 
387 IF (ITRIG.EC.1) GO TO 505 

IF (ARATIO.LT. 1.0) GC TO 460 
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IF ARATIO IS GREATER THAN ONE WHERE LIMIT = 1, THE PRESENT 

POINT VERY LIKELY IS BETWEEN TWO CLOSE ROOTS. THE FOLLOWING 

- — STATEMENTS STEP IN A SPIRAL TO LOWER ARATIO. IF ARATIO IS 

GREATER THAN ONE WHERE LIMIT = 1, POSITIVE ANALYSIS IS NOT 

— - POSSIBLE. IN AN ATTEMPT TO ANALYZE WITHOUT THE LIMIT SOME STEPS 

ARE TAKEN IN THIS CASE TOO. HOWEVER, IT IS USUALLY NECESSARY TO 

- — DIVICE CL T CNE ROOT AT THE BEST COORDINATES WHEN LIMIT = 1. 

38 5 IF (ARATIO.LT. 0.45) GO TO 490 

39 C SR AT I C = SCRT{RATI0(ND)*RATI0(ND-1) ) 

IF (SRATIO.LT » 1 . 0 T SR AT 1 0 = 1.0 
355 IF (LTRIG. EC. 1) GO TO 405 
IF (LTRIG. EQ. 2) GO TC 410 
IT (LTRIG.GT.2) GO TC 415 
XX * X 
YY * Y 

4CC LTRIG = 1 
ISIGN = C 
CAX = ARATIO 
X = X + C.l/SRATIO 
GC TC 30C 
4C5 LTRIG = 2 

CAY = ARATIO 
X = X - C.l/SRATIO 

Y = Y + C.l/SRATIO 
GC TC 30C 

41C Y = Y - C.l/SRATIO 

XCCN = (CAX - 0AY1/0AX 
YCCN (CAX - AR AT 1 0 ) /O AX 
ZCCN =; SC RT ( XCON **2 + YC0N**2) 

THETAO = ATAN2(YC0N,XC0N) - CTHETA 
XPCLE = X 
YPCLE = > 

411 THETA = CTHETA 

IF (ZCON.GT. 1.414) SRATIO = 10.0*SRATI0 
RAC * 0.6*ABS(THETA)/$RATI0 
DELTX = P AD*COS( THETA + THETAO) 

DELTY = R AC*S INC THET A + THETAO) 

412 CARAT = CAX 

413 XL A ST = X 
YLAST = V 

X = X POLE + DELTX 

Y = V POL E + DELTY 
LTRIG = LTRIG + 1 
GC TC 30C 

415 IF ( ARATIO. LE. 0.45) GO TO 490 
IF ( ARATIO. GT.OARAT) GO TO 420 
IF (LTRIG. NE. 3. OR. ARATIO. GE. 1.0) GO TO 418 
IF (SRATIO.LT. 10 00.0. OR .ARATIO.LT. OAR AT- 0*01) GO TO 418 
LTRIG = C 
GC TC 48C 

418 IF (LTRIG. GT. 101 GO TO 425 
THETA = THETA + CTHETA 
RAC - 0 • € * A BS ( THETA ) /SR AT 10 
DELTX = PAC*COS( THETA + THETAO) 

DELTY = PAD*SIN( THETA ♦ THETAO) 

CARAT = ARATIO 
GC TC 412 
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42C TF (LIMIT. GT.O.OR. LTRIG. EQ. 3 ) GO TO 425 
X = XLAST 

Y = YLAST 
DTHET A = -CTHETA 
LTRIG = -1 

GC TC 30C 
425 X = XX 

Y * YY 

IF (LIMIT. GT.O) GO TO 512 
LTRIG = 2 
ISIGN = ISIGN + 1 
IF (ISIGN - 21 430,432,434 
43C THETAO = THETAO - 3.0+DTHETA 
GC TC 411 

432 THETAC = THETAO + 6.0*DTHET A 
GC TC 411 

434 IF (ISIGN. LE. 3) GO TC 438 

THETAO = THETAO - 6.0*DTHET A 
SR AT 10 = 10.0*SftATIO 
GC TC 411 

438 THETAC = THETAO + 3.0*DTHET A 
GC TC 411 
44 C NO = ND - 1 
LIMIT = 1 

IF (JTRIG.EO.l) GO TO 450 

IF ( ABS(XP(ND-1) ).LT.0.6*T0LIQ.AND.ABS(YP(ND-1) > .LT .0. 6*T0LI Q) 

X GC TC 450 
LIMIT = C 

D = ( XP ( N D-l ) ) ** 2 + (YP(ND-1 ) >**2 

RRATIC(N'C-l) * ( XP ( NO *XP ( ND-1 )+YP(ND)*YP( ND-1) ) /0*CA ( NO ) /CA ( ND-1 ) 
CRATIC(NC-l) = ( YP( NC) *XP( ND-1) -XP(ND)*YP (ND-1) ) /D*C A ( ND ) /CA ( ND-1 ) 
D3 = (RRATIO(ND-l) )**2 + ( CRATI0(ND-1 ) )**2 
X = X - PRATIO(NC— 1)/D3 

Y * Y + CRATI0(NC-1)/D3 
45C JTPIG = 2 

GC TC 30C 
455 NC * ND - 1 
JTRTG = 2 
GC TC 39C 

46C IF (ARATIC - 0.45) 490,490,470 
47 C IF (LIMIT. EC. 0) GO TO 480 
GC TC 39C 

48C IF (ND.GT.NDEG/2 ) GO TO 498 
NC = ND + 1 
482 M = M - 1 

DC 484 1*1, M 
FK * N - NCRV - I + 1 
484 B ( I) * e(I)*FK 
GC TC 32C 

49C IF (RATIC(NC-1).LT.5.0/T0LI) GO TO 495 
X = X + C.2*T0LI 

Y = Y + C . 2*T0L I 
GC TC 3CC 

495 LTRIG = C 

IF (RATIC (NC-D.GT.l.O/TOLI ) GO TO 520 
IF (LIMIT. GT.O) GC TC 500 
498 X = X -RPATI0(ND-1)/C3 
Y= Y + CPATI0(ND-1)/C3 
GC TC 3CC 
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! 

I 50C IF < ABS (PATIC(ND-l) ) .LT •1000*01 GO TO 512 

1 XCLC = X 

? YCLC - Y 

OROTIC * RATIO(NC-l) 

ITRIG = 1 

DELTY = C.l/ORATIO 
IF (Y.LT.O.O) DELTY = -DELTY 

Y = Y + CELTY 
GO TC 3CC 

505 IF (ARATIO.GE.O.45 ) GO TO 510 
IF (LIMIT. EC. 0) GO TO 520 
IF ( ARATIO.GT.SQRT(ORATIO) ) GO TO 510 
CELTY = 2.0+DELTY 

Y = Y + CELTY 
GO TC 3QC 

51C X = XCLD 

Y = YGLD 

C CHECK RACK TO SEE IF ANOTHER REAL ROOT HAS BEEN FOUND. 

512 IF < ABS(Y) .GT.10.0*T0LI ) GO TO 515 
NC = 1 
IREAL = 2 
RETURN 

DIVIDE OUT CNLY ONE ROOT SINCE THE ROOT POINT COULD NOT BE 

LOCATED AS WELL AS DESIRED. 

515 MULT = 1 
XI = X 
YI = Y 
GC TC 54C 

52 C IF <PGLY(NC-l).LT.POLY(ND+in GO TO 525 
IF (ND.GT.2) GO TO 523 
NC = N D ♦ 1 
LTRIG = -1 
GC TC 482 
522 ND = ND - 1 
GC TC 3CC 

525 IF ( A8S(Y ) .GT .10 .0*T0L I ) GO TO 530 
IREAL = C 
XX= X 
NNC=» ND 
RETURN 
528 X = >X 
NC = NNC 
IREAL = C 

THE MULTIPLICITY OF THE ROOT AND ITS VALUE WITHIN TOLI ARF 

NCfc AT HANC. ESTIMATE THE Z LOCATION OF THE NEXT CLOSEST ROOT 

THAT IS NOT THE COMPLIMENT OF THE PRESENT ROOT. 

53C MU L T = NC- 1 
XMLL T = MILT 

IF ( NCEG - 2 *MUL T . LT • 3 ) GO TO 535 
RAT I C R= FRATIO(NC)/(XMULT + 1.0) 

R A T I C 1= CRATIG(NCf/( XMULT + 1.0) + XMULT/I 2.0*Y) 

OX => RATIOR*Y/(Y *( RAT I OR** 2 + RATI0I**2) + RATIOI) 
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XI 3 X - DX 

YI » SORT (ABS((RATIOR*Y + RATIOI*OX)*< DX*DX+ Y*Y ) / ( RAT I OR*Y - 
X RATIOI*CX) ) ) 

535 D2 = ( RRAT I C ( MULT ) ) **2 + (CRATIOl MULT ) )**2 
X= X - PRATIO(MULT)/C2 
Y= Y + CR AT 10 (MULT ) /02 

— NUMEER THE ROOTS AND DIVIDE THEM OUT OF THE POLYNOMIAL FIRST 

AS COMPLEX NUMBERS TO FIND THE COMPLEX RESIDUAL. THEN GO BACK AND 

--- DIVIDE CLT THE ROOT AND ITS COMPLIMENT TOGETHER SO THAT ONLY RFAL 

NUMBER ALGEBRA IS USED. 

54C DC 565 K*1,MULT 
DC 545 1*1 *N 
RA ( T ) = A(I) 

545 CA ( I ) = C .0 
J + l 

XRCCTU) = X*SCALE 
YRCCT(J) =» Y*SCALE 
N= N - 1 

DC 550 I *2 * N 

XPN* R A ( I ) + X*R A ( 1-1 ) - Y*C A ( I- 1 ) 

CA(d= C A ( I ) + X*CA ( I— 1 1 + Y*R A ( I- 1 ) 

55C RA ( I ) n XRN 

XR(J) = ( R A ( N + l ) + X*RA(N) - Y*CA( N ) ) /A (N+l ) 

555 YR ( J ) = (CACN+1) * X*CA(N) + Y*RA ( N ) ) / A ( N + l ) 

J= J + l 

XRCCT(J) * X*SCALE 
YRCCTU) *-Y*SCALE 
N = N - 1 

A ( 2 ) = A ( 2 ) + 2.0 *A ( 1 ) *X 
IF (N.EC.l) GO TO 562 
K= N + I 

DC 560 1*3, K 

56C A(I)* A(I) ♦ 2.0 *A (I- 1 ) *X - A(I-2)*(X**2 + Y**2) 

562 XR ( J ) = 1.0 + ( A ( N+l ) *X - A(N)*(X**2 + Y**2 ) ) /A ( N+2 ) 

565 YR(J)n 0.0 
X* XI 
Y* YI 

IF (X.GT.10.0) X = 10.0 

IF ( (N+1J/2.E0.N/2) IREAL = 1 

RETURN 

ENC 


SUERCUT IN E CPOLY 

C THIS IS THE BASIC ROUTINE FOR FINDING THE VALUE OF THE 

C POLYNOMIAL AND ITS DERIVATIVES IN COMPLEX NUMBERS. 

DCLELE PRECISION A, ARATIO, B, CA, CRATIO, D, D2, FK, P, POLY, Q, 
X RA, RATIO, RRAT 10 , SCALE, X, XL AST, XP, XPN, XPOLY, XR, 

X XRCCT, V, YP, Y POLY , YR, YRCOT 

CCMMCN /FOLYN/ ARATIO, D, DELTX , DX, 02, FK, I, INK, IREAL, 

1 ITRIG, IVI, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N, NC , ND, 

2 NCEG , NCRV, NEXP, NEXPS, NK, NM, NUM, OARAT, Q, SCALE, SRATIO, 
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3 TCLI, TCLIQ, TOLR, TOLRQ, X, XI, XL AST , XMULT » XOLO, XPN, Y, 

4 A ( 1 CO ) , B ( 100) , CA(IOO), CRATIO(IOO), P(100), POLY(IOO), RA(IOO), 

5 RATIC(ICO), RRATIOC 100) , XP(IOO), XPOLY(IOO), XR(IOO), XROOT(IOO) 

6 , YP ( 100 ) » YPOLY(IOO), YR(IOO), YROOT(IOO) 

32 C XP(NCRV+11= 1.0 

YP ( N DRV + 1 ) 3 0.0 
L= P 

IF 1P.EC.0) GO TO 345 
NC = NEXPS - NEXP 
IF (L.EC.N-1) GO TO 328 
NF 3 NC 

325 IF (ABS(B(L+1)/CA(NDRV> ) .GT • 0.00001*8. 0**NM*T0L IQ) GO TO 328 
L- L-l 

IF (L.EC.O) GO TO 340 
NM NM 4 NC 
GC TC 325 
32E I = 1 

33C IF (I. EC. L+l) GO TO 340 

XP (NCRV+1 ) 3 XP(NDRV+l)*Bm 
II 3 I 
NP = NC 

332 IF (ABSIBC I+1)/B(II) ).GT. 0.00001*8. 0**NM*T0LIQ) GO TO 335 
XPN 3 XP ( N DRV+ 1 ) *X - YP( NDRV+1 ) *Y 
YP(NCRV+1> 3 XP(NCRV+1)*Y + YP(NDRV+1)*X 
XP (NCRV+1 )= XPN 
IF < I .EC.L ) GO TO 338 
I 3 I + 1 
NM NM 4 NC 
GC TC 332 

335 XPN 3 XP(NDRV+1)*X - YP(NDRV+1)*Y 

YP ( NCRV + 1 ) = XPl NDRV+1 ) *Y + YP(NDRV+1)*X 
XP ( NCRV+1 ) 3 XPN/B ( I + 1 ) + 1.0 
1=1 + 1 
GC TC 33C 

338 XP(NCRV+1) 3 XP(NDRV+1)/B( I +1 ) + 1.0 
34 C IF (L.EC.P) GO TC 345 

YPCNCRV+l) 3 YP( NCRV+1) / R ( L+ 1 ) 

PL = P - L 
DC 342 I 3 1 , PL 

XPN = XP(NDRV+1)*X - YP ( NDRV+1 ) * Y 
YP (NCRV + 1 ) = XPt NDRV + 1 ) *Y + YP(NDRV+1)*X 
342 XPtNCRV+1) = XPN 

YPCLYCNCFV+l) = YP l NDRV + 1 ) *B ( L + l ) 

GC TC 348 

345 YPCLYfNCFV+1) = YPCNDRV+1) 

YP ( NC RV + 1 ) = YP( NDRV+1 )/B( L + l ) 

348 XPCLY (NCPV + 1) = XP [ NDRV + 1 ) *B ( L + 1 ) 

CA ( N CRV+ 1 ) = B( L + 1 ) 

NCR V= NDFV+1 
IF (NCRV.GT.ND) RETURN 
348 P= P-1 

DC 350 1*1, P 
FK= N-NCRV-I + 1 
35 C BID* B { I >*FK 
GC TC 32C 
ENC 
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SUBROUTINE SCALER 

THIS SUBROUTINE SCALES THE POLYNOMIAL SO THE LOCAL ROOT IS 
OF CRCER ONE. 

DOUBLE PRECISION A, ARATIO, B, CA, CRATIO, 0, D2, FK, P, POLY, Q, 

X RA, RATIO, RRATIO, SCALE, X, XL AST , XP, XPN, XPOLY, XR, XROOT, Y, 
X YF, YPCLY, YR, YROOT 

CCMMCN /FOLYN/ ARATIO, D, DELTX, DX, D2, FK, I, INK, IREAL, 

1 ITRIG, IW, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N, NC , NO, 

2 NCEG , NCRV, NEXP, NEXPS, NK , NM, NUM, OARAT, Q, SCALE, SRATIO, 

3 TCLI, TCLIG, TOLR , TOLRQ, X, XI, XL AST , XMULT , XOLO, XPN, Y, 

4 A (ICC ) , B ( 100) , CA(IOO), CRATIO(IOO), P(100), POLY(IOO), RAC100), 

5 RATIC(ICO), RRATIO(IOO), XP(IOO), XPOLY(IOO), XR(ICO), XROOT(IOO) 

6 ,YP ( ICO ) , YPOLY ( 100 ) , YR(IOO), YROOT(IOO) 

NK * C 

INK = 1 

IF (IREAL. LE.O) GO TC ’.0 
XPN = 0.36067 376*AL0G(ABS(X) ) 

GC TC 20 

1C XPN =i 0.18033 688*AL0G(X**2 + Y**2) 

2C IF (XPN. GT. 0.0) GO TO 30 
XPN = -XFN 
INK * -1 
3C DC 4C 1=1,6 

TF (XPN.LT.0.5.0R.IABS(NC + (NDEG* ( NK + INK ) ) /2 ) .GT . 20 ) GO TO 50 
NK = NK 4 INK 
4C XPN = XPN - 1.0 
5C IF (NK. EC .0 ) RETURN 
NC = (N0EG*NK)/2 
NEXP = NEXP + NK 
SCALE = 16.0DO**NEXP 
DC 60 1 = 1, N 
NM = NK*(I-1) - NC 
6C AID ^ A ( I ) /16.0C0**NM 
RETLRN 
ENC 


SUB R CUT I N E OED 

C IN THIS SUBROUTINE THE LAST TWO ROOTS ARE CALCULATED DIRECTLY 

C ANC ALL THE ROOTS ARE PRINTED. 

DCIELE PRECISION A, ARATIO, B, CA, CRATIO, 0, D2, FK, P, POLY, Q, 

X RA, RATIO, RRATIO, SCALE, X, XL AST, XP, XPN, XPOLY, XR, 

X XRCCT , Y , YP, YPOLY, YR, YROOT 

CCMMCN /FOLYN/ ARATIO, D, DELTX, DX, D2, FK, I, INK, IREAL, 

1 ITRIG, IW, J, JTRIG, L, LIMIT, LTRIG, M, ML, MULT, N, NC , ND, 

2 NCEG, NCRV, NEXP, NEXPS, NK , NM, NUM, OARAT, Q, SCALE, SRATIO, 

3 TCLI, TCLIC, TOLR, TOLRQ, X, XI, XL AST, XMULT, XOLD, XPN, Y, 

4 A ( ICC ) , B(100), CA(IOO), CRATIO(IOO), P(100), POLY(IOO), RAUOO), 

5 RATIC(ICO), RRATIOUOO), XP(IOO), XPOLYUOO), XR(IOO), XROOTUOO) 
6 , YP ( 100 ) , YPOLY ( 100 ) , YR(IOO), YROOT(IOO) 

57C IF (N-2) 600,575,580 
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n o 


575 X = -A(2)/A(l) 

J = J + 1 

XRCCT(J) a X+SCALE 
YRCCTtJ) = 0.0 
XR ( J ) = C.O 
YR ( J J = C.O 
GC TC 60C 

WHEN THERE ARE ONLY TWO ROOTS LEFT, IT IS EASIER TO SOLVE FOR 

THEM DIRECTLY. 

58C Q = A(2)«*2 - 4. 0* A ( 1 ) * A ( 3 ) 

IF (C.LT.O.O) GO TO 590 
Q * CSQRTtC) 

J * J + 1 

XRCCT(J) » (Q - At 2 ) )/( 2.0*A(1) )*SCALE 
YRCCT(J) a 0.0 
XR ( J I = C.O 
YR ( J ) = C.O 
J = J + I 

XRCCTtJ) a -{Q +■ A<2))/(2.0*Atl))*SCALE 
YRCCTtJ) * 0.0 
XR ( J ) = C.O 
YR ( J ) = C.O 
GC TC 6CC 
59C J a J + 1 

XRCCT(J) a -A(2) /(2.0*A( 1) )*SCALE 
YRCCTtJ) * DSQRT (-Q ) / ( 2.0*A( 1) )*SCALE 
XR ( J ) = C.O 
YRtJ) = C.O 
J = J + 1 

XRCCTtJ) a XROOTtJ-1) 

YRCCTtJ) = -YROOTIJ-1) 

XR C J ) a C.O 
YRtJ) = C.O 
60C WRITE t IV » 1040) 

DC 610 J»1,NUM 

61C WRITE (IV, 1050) J, XROOT(J), YROOT(J), XR ( J ) , YRtJ) 

WRITE (IV, 1060) 

104C FORMAT (////// 39X , 54( 1 H* ) / 39X , 1H*, 52X, 1H*/39X , 1H* , 4X, 48H*** THE 

1 CCMPUT EC ROOTS OF THE POLYNOMIAL *** , 1H* /39X, 1H* , 52X , 1H* /I 3X , 

2 106 t 1H*)/13X».1H*»8X»1H*,47X»1H*»47X*1H*/13X»1H*,8X,1H**47X,1H*, 

3 11X,25HF0LYN0MI AL REMAINDER WHEN, 11X, 1H*/13X, 10H* ROOT *,8X, 

4 9HREAL F ART , 13X , 9H I M AG INARY, 8X , 1H* , 13X , 21HR00T WAS DIVIDED OUT , 

5 13X,1H*/13X» 1H* , IX, 7HN UMBER ,1H*,9X,7H0F ROOT, 13X, 12HPART OF ROOT 

6 »6X»1H*,47X,1H*/13X»1H*,8X,1H*,47X,1H*,8X, 9HRE AL P ART , 1 1 X , 14H I MAG 

7 I NARY PAPT,5X,1H*/13X,1H*,8X,1H*,47X, 1H*, 47X, 1H* /13 X , 106 ( 1H* ) /13X , 

8 IF*»8X,1H*, 47X, 1H* » 47X , 1H* ) 

105C FORMAT (13X,1H*, 2X , 13 , 3X , 1H* , 2022. 1 2 , 3X , 1H* , 2022 . 12 , 3X , 1H* ) 

106C FORMAT ( 1 3X , 1H*, 8X, 1H*, 47X , 1H*, 47X, 1H*/1 3X , 106( 1H* ) ) 

RETLRN 

ENC 
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CO 


Examples 


THE INPUT POLYNOMIAL COEFFICIENTS, THEY ARE IN ORDER IF READ IN ROWS. 

O.IOOOOCOOCCOOOOOOD 01 -0. 49998999999999990 02 0. 9999600000000000D 03 -0. 99993999999999990 04 0.4999600000000000D 05 

-C.9999CCOGCCOOOOOOD 05 


************************** ************* ******** ******* 

* * 

* *** THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 

** 4 ******************************************************************************************************* 


* 


* 




* 



* 

* 


* 




* 

POLYNOMIAL 1 

REMAINDER WHEN 

* 

* 

RCOT 

* 

REAL PART 


IMAGINARY 

* 

ROOT WAS ' 

DIVIDEO OUT 

* 

* 

AUMBER 

* 

OF ROOT 


PART OF ROOT 

* 



* 

* 


* 




* 

REAL PART 

IMAGINARY PART 

* 

* 


* 




* 



* 

************ 

********************************************************************************************** 

* 


* 




* 



* 

* 

1 

* 

0.100000C000Q5D 

02 

0. 

* 

0. 8881784197000-15 

0. 

* 

* 

2 

* 

0.100000C00005D 

02 

0. 

* 

-0.1221245327090-14 

0. 

* 

* 

3 

* 

0 .1000000000050 

02 

0. 

* 

-0.9992007221630-15 

0. 

* 

* 

4 

♦ 

0 .1000000000050 

02 

0. 

* 

0.2053912595560-13 

0. 

* 

* 

5 

* 

0.9998999998 120 

01 

0. 

* 

0. 

0. 


* 


♦ 




* 



* 


*************** ****4 ************************************************************** ********************* 


4 

Exact polynomial: P(x) = (x - 10) (x - 9. 999) 

Comment: This problem illustrates the point that good resolution capability is maintained near multiple roots when the multi- 
ple root in a group is found first. 



THE INPUT POLYNOMIAL COEFFICIENTS 


THEY ARE IN ORDER IF READ IN ROWS 


O.IOOCCCOOOCOOOOOOD 01 -0.4999990000000000D 01 0.89999600000000000 01 -0. 6999950000000000D 01 0. 1999980000000000D 01 


****************************************************** 

* * 

* *** THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 
********************************************************************************************************** 


* 


* 





* 



* 

* 


* 





* 

POLYNOMIAL 

REMAINDER WHEN 

* 

* 

RC0T 

* 

REAL PART 



IMAGINARY 

* 

ROOT WAS 

DIVIDED OUT 

* 

* 

AUMBER 

* 

OF ROOT 



PART OF ROOT 

* 



* 

* 


* 





+ 

REAL PART 

IMAGINARY PART 

* 

* 


* 





* 



* 

********************************************************************************************************** 



* 





* 



♦ 


1 

t 

0.999999999980D 

00 

0. 


* 

-0. 

0. 

* 


2 

* 

0 .9999999999800 

00 

0. 


* 

-0.6661338147750-15 

0. 

* 


3 

* 

0 .2000000000000 

01 

0. 


* 

0. 

0. 

* 


4 

* 

0.9999900000390 

00 

0. 


* 

0. 

0. 

* 



♦ 





* 



* 


********************************************************************************************************** 


Exact polynomial: P(x) = (x - l) 2 (x - 0. 99999) (x - 2) 

Comment: This problem again shows good resolution near a multiple root. However, in this case the decision between 
x = 1.000000 and x = 0. 999993 for the location of the multiple root was barely possible or maybe even lucky. 
When the correct choice was made, the computed root values were quite good. 



THE INPUT POLYNOMIAL COEFFICIENTS. THEY ARE IN ORDER IF READ IN ROWS 


O.1OOCOCOOOCOOOO0OD 01 “0. 50000C9999999999D 01 


0.9000040000000000D 01 -0. 7000050000000000D 01 0. 200D020000000000D 01 


** ***** ************ ******** ******** ******* ************ 

* * 

* **+ THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 

********** ************************************************************************************************ 


* 


* 




* 



* 

* 


* 




* 

POLYNOMIAL 

REMAINDER WHEN 

* 

* 

RCCT 

* 

REAL PART 


IMAGINARY 

* 

ROOT WAS 

DIVIDED OUT 

* 

* 

rUMRER 

* 

OF ROOT 


PART OF ROOT 

* 



* 

* 


* 




* 

REAL PART 

IMAGINARY PART 

* 

* 


• 




* 



* 

* 44 * ****** ****** ******* ****** * ********************************************************************* ******* 

* 


* 




* 




* 

1 

* 

0.100000666680D 

01 

0. 

* 

-0.3 33066907 3880- 15 

0 . 


* 

2 

* 

0. 100 000 66668 0D 

01 

0. 

* 

0.18873 79 14186 D- 14 

0. 


* 

3 

* 

O. 2 OOOOOCOOOOOD 

01 

0. 

» 

0 . 

0. 


* 

4 

* 

0.999996666400D 

00 

0. 

* 

0. 

0. 


* 


* 




* 




** 4 ************************************************************************************ ******************* 


Exact polynomial: P(x) = (x - l) 2 (x - 1.00001)(x - 2) 

Comment: This problem is similar to the previous one* In this case the wrong choice between x = 1.000000 and 

x = 1.0000067 was made for the location of the multiple root. Of course, the root values are not as good as 
they would have been if the correct decision could have been made; but note that the centroid of the three roots 
in the group is quite accurate. 



THE INPUT POLYNOMIAL COEFFICIENTS 


THEY ARE IN OROER IF READ IN ROWS 


O.IOCOOCOOQCCOOOOOD 01 -0. AOOOOCOOOOOOOOOOD 01 0. 60CC001999999999D 01 -0. AOOOOOAOOOOOOOOOD 01 0* 100D00200000 1000D 01 


****** ** *************************** ********** ********* 


* 


* 


* *** THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 

*4 ********************************************************************************** ********************* 


* 


* 


* 


* 


* RCOT * 

* MJMBER * 

* * 


REAL PART 
OF ROOT 


* 

IMAGINARY * 

PART OF ROOT * 

* 


POLYNOMIAL REMAINDER WHEN * 

ROOT WAS DIVIDED OUT * 

* 

REAL PART IMAGINARY PART * 


% * * * 
********************************************************************************************************** 


1 

* 

* 

0*999999999 8 AO D 

00 

0*9999999993330-03 

* 

* 

-0*9002862256570-16 

0.A3368D00163AD-18 

* 

* 

2 

* 

0.9999999998AOD 

00 

-0* 9999999 9933 3D-03 

* 

-0.111022302A63D-15 

0. 

* 

3 

* 

0*9999 99 9998 AOD 

00 

0 • 999999999333D-03 

* 

0.3 1390 38 91732 D-1A 

-0.6A0736713725D-12 

* 

A 

* 

* 

0*9999999998 A00 

00 

-0*9999999993330-03 

* 

* 

0* 32196A6771A1D-1 A 

0* 

* 

* 


********************************************************************************************************* 


Exact polynomial: P(x) = (x - 1 + 0.001i)^(x- 1 - O.OOli)^ 

Comment: This problem is an illustration of multiple roots in the complex plane. The problem also shows capability to dis- 
tinguish complex roots from nearly real roots. 



-a 

05 


THE INPUT POLYNOMIAL COEFFICIENTS, THEY ARE IN ORDER IF READ IN ROWS. 

G.IOOCCCCOOCOCOOOOD 01 -0. 8000000000000000D 01 0. 2 500000199999999D 02 -0. 4000001200000000D 02 0 . 350D002400000100D 02 

-C • I60CC0 20CC 000399D 02 0. 30000C60000029990 01 


****************************************************** 

* * 

* *** THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 
********************************************************************************************************** 


• 


* 




* 



* 

* 


♦ 




* 

POLYNOMIAL ! 

REMAINDER WHEN 

* 

* 

RCCT 

* 

REAL PART 


IMAGINARY 

* 

ROOT WAS 1 

DIVIDED OUT 

* 

* 

hUMBER 

* 

OF ROOT 


PART OF ROOT 

* 



* 

* 


* 




* 

REAL PART 

IMAGINARY PART 

* 

* 


* 




* 



* 

** i ******************************************************************************************************* 

* 


* 




* 



* 

* 

1 

* 

0.9990765817 760 

00 

0. 

* 

0.2442490654180-14 

0. 

* 

* 

2 

* 

0 .999687 7264080 

00 

0.1487630401670-02 

* 

-0.3314258028320-13 

0.713466965356 D- 16 

* 

* 

3 

* 

0.9996877264080 

00 

-0.1487630401670-02 

* 

-0.3430589146090-13 

0. 

* 

* 

4 

* 

C.3O0OOOCOOO0OO 

01 

0. 

* 

-0.1998401444330-14 

0. 

* 

* 

5 

* 

0.1000773982700 

01 

0.9538712732370-03 

* 

0. 

0. 

* 

♦ 

6 

* 

0.1000773982700 

01 

-0.9538712732370-03 

* 

0. 

0. 

* 

* 


* 




* 



* 


** 4 ******************************************************************************************************* 


Exact polynomial: P(x) = (x - l)(x -1 + 0. 001 i) 2 (x - 1 - 0. 001 i)^(x - 3) 

Comment: This is the previous problem with a real root between the complex roots. If the program could get the multiple 
root first, the resolution would be quite good; but the program is structured to find real roots first. Since the 
polynomial was found to change sign in the vicinity of the root group near x = 1. 0, a real root is known to exist. 
The real root could not be resolved with the ratio of derivatives theory so a root was taken out near the center of 
the group. 



THE INPUT POLYNOMIAL COEFFICIENTS 


THEY ARE IN ORDER IF READ IN ROWS 


O.IOCCOCOOOCOOOOOOO 01 -0. -0.100C001000001000D 07 -0. 0.1003001000001000D 07 

-C. -O.IOOOOCOOOOOOOOOOD 01 


POLYNOMIAL REMAINDER WHEN 
ROOT WAS DIVIDED OUT 


****************************************************** 

* * 

* *** THE COMPUTED ROOTS OF THE POLYNOMIAL *** * 

* * 

** 4 ****************************************************************************************************** 
* * * 

* * * 

* RCOT * REAL PART IMAGINARY * 

* NUMBER * OF ROOT PART OF ROOT * 

* * * 

* * * 

*44 ******* **** 4 v** *************************************************************** ************************ 

* 

0 . * 

0 . * 

0 . * 

0 . * 

0 . * 

0 . * 


REAL PART 


IMAGINARY PART 


* 

* 

* 

* 

* 

* 

* 

* • • 

** 4 41****** *************************************************** ******************************************** 


O.lOOOOOCOOOOOD-02 
-O.lOOOOOCOOOOOD-02 
-O.IOOOOOCOOOOOD 01 
0.100000000000D 01 
0.100000000000D OA 
-O.IOOOOOCOOOOOD OA 


-0.781597009336D-13 
-0.2220AA60A925D-15 
-0. 2220AA60A925D-1 A 
-0. 1 1 102 2302 A63D-15 

0. 

0 . 


0. 

0. 

0. 

0. 

0 . 

0 . 


Exact polynomial: P(x) = (x + 0. 001) (x - 0. 001)(x + l)(x - l)(x + 1000)(x - 1000) 

Comment: This problem shows good root accuracy over a wide range of root size. This is accomplished by (1) polynomial 
scaling and (2) extracting the lower magnitude roots first. 



APPENDIX D 


OUTLINE OF THE BASIC LOGICAL STEP PROCESS FOR THE RATIO 
OF DERIVATES ANALYSIS FOR REAL ROOTS 

(1) First x try at a given root 

ND = 2 

If |P(1)| < TOLRQ and |P(2)| < 10 x TOLR, ND = ND + 1. 

(2) Each new x try at a given root 

(a) ARATIO < 0. 45 (Appear to be at right root multiple, (ND - 1), and zeroing in. ) 

(1) If | RATIO(ND - 1) | > 5/TOLR, too close to the root for good analysis. 

Back away to get in the 1/TOLR to 5/TOLR range. 

(2) If 1/TOLR < |RATIO(ND - 1) | < 5/TOLR, make the final root correction 

and divide the root out of the polynomial. 

(3) If 1/TOLR > |RATIO(ND - 1) | 

(a) LIMIT = 0, ( | P(ND - 1) | is valid) make x correction by 

l/RATIO(ND - 1) to get closer to the root. 

(b) LIMIT = 1, ( |P(ND - 1) | is not valid. It is set to TOLRQ, and it is 

assumed to have greater magnitude than the actual value. ) Set 
KTRIG = 1. 

(1) If | RATIO(ND - 1) | < 1000, use present x value as root 

(2) If |RATIO(ND - 1) | > 1000, step away from the present x to a 

point where a valid RATIO(ND - 1) is reached so the 
l/RATIO(ND - 1) final correction can be made to give a better 
root value. 

(b) 0. 45 < ARATIO < 1. 0. (Either at too low ND or too far from the root yet. ) 

( 1) If ND < NDEG. Increase ND by 1 and come through the analysis again at 

the same x. 

(2) If ND = NDEG. Not yet near enough to root. Step closer by l/RATIO(ND) 

correction. Set JTRIG = 2. 

(c) ARATIO > 1. 0. (In general, between roots. ) From figures 2 and 3, it can be 

seen that root candidates occur on either side of such an x location. Set 
IMX = 1 to indicate the first of the two is being investigated by stepping along 
the real axis with LTRIG the step counter. Step with doubled step size until: 
(1) ARATIO < 0. 45 

(a) LIMIT = 0. This point is a root. Do normal analysis. 

(b) LIMIT = 1. Check the other point of symmetry in the effort to try to 

determine which is the root. 
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(2) LTRIG = 8 

(a) IMSURE = 0. Give up the search for a real root here. 

(b) IMSURE = 1. A real root is known to exist in this region. Continue 

stepping until a XHIGH or XLOW limit is reached and then step in 
the other direction. 
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